diff options
Diffstat (limited to 'ccx_2.14_linpack.f')
| -rw-r--r-- | ccx_2.14_linpack.f | 218 |
1 files changed, 218 insertions, 0 deletions
diff --git a/ccx_2.14_linpack.f b/ccx_2.14_linpack.f new file mode 100644 index 000000000000..c737987114ad --- /dev/null +++ b/ccx_2.14_linpack.f @@ -0,0 +1,218 @@ +c +c L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License” +c or “3-clause license”) +c Please read attached file License.txt +c + subroutine dpofa(a,lda,n,info) + integer lda,n,info + double precision a(lda,*) +c +c dpofa factors a double precision symmetric positive definite +c matrix. +c +c dpofa is usually called by dpoco, but it can be called +c directly with a saving in time if rcond is not needed. +c (time for dpoco) = (1 + 18/n)*(time for dpofa) . +c +c on entry +c +c a double precision(lda, n) +c the symmetric matrix to be factored. only the +c diagonal and upper triangle are used. +c +c lda integer +c the leading dimension of the array a . +c +c n integer +c the order of the matrix a . +c +c on return +c +c a an upper triangular matrix r so that a = trans(r)*r +c where trans(r) is the transpose. +c the strict lower triangle is unaltered. +c if info .ne. 0 , the factorization is not complete. +c +c info integer +c = 0 for normal return. +c = k signals an error condition. the leading minor +c of order k is not positive definite. +c +c linpack. this version dated 08/14/78 . +c cleve moler, university of new mexico, argonne national lab. +c +c subroutines and functions +c +c blas ddot +c fortran sqrt +c +c internal variables +c + double precision ddot,t + double precision s + integer j,jm1,k +c begin block with ...exits to 40 +c +c + do 30 j = 1, n + info = j + s = 0.0d0 + jm1 = j - 1 + if (jm1 .lt. 1) go to 20 + do 10 k = 1, jm1 + t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1) + t = t/a(k,k) + a(k,j) = t + s = s + t*t + 10 continue + 20 continue + s = a(j,j) - s +c ......exit + if (s .le. 0.0d0) go to 40 + a(j,j) = sqrt(s) + 30 continue + info = 0 + 40 continue + return + end + +c====================== The end of dpofa =============================== + + subroutine dtrsl(t,ldt,n,b,job,info) + integer ldt,n,job,info + double precision t(ldt,*),b(*) +c +c +c dtrsl solves systems of the form +c +c t * x = b +c or +c trans(t) * x = b +c +c where t is a triangular matrix of order n. here trans(t) +c denotes the transpose of the matrix t. +c +c on entry +c +c t double precision(ldt,n) +c t contains the matrix of the system. the zero +c elements of the matrix are not referenced, and +c the corresponding elements of the array can be +c used to store other information. +c +c ldt integer +c ldt is the leading dimension of the array t. +c +c n integer +c n is the order of the system. +c +c b double precision(n). +c b contains the right hand side of the system. +c +c job integer +c job specifies what kind of system is to be solved. +c if job is +c +c 00 solve t*x=b, t lower triangular, +c 01 solve t*x=b, t upper triangular, +c 10 solve trans(t)*x=b, t lower triangular, +c 11 solve trans(t)*x=b, t upper triangular. +c +c on return +c +c b b contains the solution, if info .eq. 0. +c otherwise b is unaltered. +c +c info integer +c info contains zero if the system is nonsingular. +c otherwise info contains the index of +c the first zero diagonal element of t. +c +c linpack. this version dated 08/14/78 . +c g. w. stewart, university of maryland, argonne national lab. +c +c subroutines and functions +c +c blas daxpy,ddot +c fortran mod +c +c internal variables +c + double precision ddot,temp + integer case,j,jj +c +c begin block permitting ...exits to 150 +c +c check for zero diagonal elements. +c + do 10 info = 1, n +c ......exit + if (t(info,info) .eq. 0.0d0) go to 150 + 10 continue + info = 0 +c +c determine the task and go to it. +c + case = 1 + if (mod(job,10) .ne. 0) case = 2 + if (mod(job,100)/10 .ne. 0) case = case + 2 + go to (20,50,80,110), case +c +c solve t*x=b for t lower triangular +c + 20 continue + b(1) = b(1)/t(1,1) + if (n .lt. 2) go to 40 + do 30 j = 2, n + temp = -b(j-1) + call daxpy(n-j+1,temp,t(j,j-1),1,b(j),1) + b(j) = b(j)/t(j,j) + 30 continue + 40 continue + go to 140 +c +c solve t*x=b for t upper triangular. +c + 50 continue + b(n) = b(n)/t(n,n) + if (n .lt. 2) go to 70 + do 60 jj = 2, n + j = n - jj + 1 + temp = -b(j+1) + call daxpy(j,temp,t(1,j+1),1,b(1),1) + b(j) = b(j)/t(j,j) + 60 continue + 70 continue + go to 140 +c +c solve trans(t)*x=b for t lower triangular. +c + 80 continue + b(n) = b(n)/t(n,n) + if (n .lt. 2) go to 100 + do 90 jj = 2, n + j = n - jj + 1 + b(j) = b(j) - ddot(jj-1,t(j+1,j),1,b(j+1),1) + b(j) = b(j)/t(j,j) + 90 continue + 100 continue + go to 140 +c +c solve trans(t)*x=b for t upper triangular. +c + 110 continue + b(1) = b(1)/t(1,1) + if (n .lt. 2) go to 130 + do 120 j = 2, n + b(j) = b(j) - ddot(j-1,t(1,j),1,b(1),1) + b(j) = b(j)/t(j,j) + 120 continue + 130 continue + 140 continue + 150 continue + return + end + +c====================== The end of dtrsl =============================== + + |