1 ! from expokit
, free
for non
-commercial use
3 * This is a lightweight substitute
to the
external LAPACK routines
4 * used by EXPOKIT
. It is supplied
to ensure that EXPOKIT is
5 * self
-contained and can still run
if LAPACK is not yet installed
6 * in your environement
.
7 *----------------------------------------------------------------------|
8 subroutine DGESV
( N
, M
, A
,LDA
, IPIV
, B
,LDB
, IFLAG
)
9 integer N
, M
, LDA
, LDB
, IPIV
(N
), IFLAG
10 double precision A
(LDA
,N
), B
(LDB
,M
)
11 call DGEFA
( A
,LDA
, N
, IPIV
, IFLAG
)
12 if ( IFLAG
.ne
.0 ) stop "Error in DGESV (LU factorisation)"
14 call DGESL
( A
,LDA
, N
, IPIV
,B
(1,j
), 0 )
17 *----------------------------------------------------------------------|
18 subroutine DSYSV
(UPLO
, N
,M
, A
,LDA
, IPIV
, B
,LDB
, WRK
,LWRK
, IFLAG
)
20 integer N
, M
, LDA
, LDB
, LWRK
, IFLAG
, IPIV
(N
)
21 double precision A
(LDA
,N
), B
(LDB
,M
), WRK
(LWRK
)
22 call DSIFA
( A
,LDA
, N
, IPIV
, IFLAG
)
23 if ( IFLAG
.ne
.0 ) stop "Error in DSYSV (LDL' factorisation)"
25 call DSISL
( A
,LDA
, N
, IPIV
, B
(1,j
) )
28 *----------------------------------------------------------------------|
29 subroutine ZGESV
( N
, M
, A
,LDA
, IPIV
, B
,LDB
, IFLAG
)
30 integer N
, M
, LDA
, LDB
, IPIV
(N
), IFLAG
31 complex*16 A
(LDA
,N
), B
(LDB
,M
)
32 call ZGEFA
( A
,LDA
, N
, IPIV
, IFLAG
)
33 if ( IFLAG
.ne
.0 ) stop "Error in ZGESV (LU factorisation)"
35 call ZGESL
( A
,LDA
, N
, IPIV
,B
(1,j
), 0 )
38 *----------------------------------------------------------------------|
39 subroutine ZHESV
(UPLO
, N
,M
, A
,LDA
, IPIV
, B
,LDB
, WRK
,LWRK
, IFLAG
)
41 integer N
, M
, LDA
, LDB
, LWRK
, IFLAG
, IPIV
(N
)
42 complex*16 A
(LDA
,N
), B
(LDB
,M
), WRK
(LWRK
)
43 call ZHIFA
( A
,LDA
, N
, IPIV
, IFLAG
)
44 if ( IFLAG
.ne
.0 ) stop "Error in ZHESV (LDL' factorisation)"
46 call ZHISL
( A
,LDA
, N
, IPIV
,B
(1,j
) )
49 *----------------------------------------------------------------------|
50 subroutine ZSYSV
(UPLO
, N
,M
, A
,LDA
, IPIV
, B
,LDB
, WRK
,LWRK
, IFLAG
)
52 integer N
, M
, LDA
, LDB
, LWRK
, IFLAG
, IPIV
(N
)
53 complex*16 A
(LDA
,N
), B
(LDB
,M
), WRK
(LWRK
)
54 call ZSIFA
( A
,LDA
, N
, IPIV
, IFLAG
)
55 if ( IFLAG
.ne
.0 ) stop "Error in ZSYSV (LDL' factorisation)"
57 call ZSISL
( A
,LDA
, N
, IPIV
, B
(1,j
) )
60 *----------------------------------------------------------------------|
61 *----------------------------------------------------------------------|
62 *----------------------------------------------------------------------|
63 subroutine dgefa
(a
,lda
,n
,ipvt
,info
)
64 integer lda
,n
,ipvt
(n
),info
65 double precision a
(lda
,n
)
67 c dgefa factors a double precision matrix by gaussian elimination.
69 c dgefa is usually called by dgeco, but it can be called
70 c directly with a saving in time if rcond is not needed.
71 c (time for dgeco) = (1 + 9/n)*(time for dgefa) .
75 c a double precision(lda, n)
76 c the matrix to be factored.
79 c the leading dimension of the array a .
82 c the order of the matrix a .
86 c a an upper triangular matrix and the multipliers
87 c which were used to obtain it.
88 c the factorization can be written a = l*u where
89 c l is a product of permutation and unit lower
90 c triangular matrices and u is upper triangular.
93 c an integer vector of pivot indices.
97 c = k if u(k,k) .eq. 0.0 . this is not an error
98 c condition for this subroutine, but it does
99 c indicate that dgesl or dgedi will divide by zero
100 c if called. use rcond in dgeco for a reliable
101 c indication of singularity.
103 c linpack. this version dated 08/14/78 .
104 c cleve moler, university of new mexico, argonne national lab.
106 c subroutines and functions
108 c blas daxpy,dscal,idamax
113 integer idamax
,j
,k
,kp1
,l
,nm1
116 c gaussian elimination with partial pivoting
120 if (nm1
.lt
. 1) go to 70
124 c find l = pivot index
126 l
= idamax
(n
-k
+1,a
(k
,k
),1) + k
- 1
129 c zero pivot implies this column already triangularized
131 if (a
(l
,k
) .eq
. 0.0d0
) go to 40
133 c interchange if necessary
135 if (l
.eq
. k
) go to 10
141 c compute multipliers
144 call dscal
(n
-k
,t
,a
(k
+1,k
),1)
146 c row elimination with column indexing
150 if (l
.eq
. k
) go to 20
154 call daxpy
(n
-k
,t
,a
(k
+1,k
),1,a
(k
+1,j
),1)
163 if (a
(n
,n
) .eq
. 0.0d0
) info
= n
166 *----------------------------------------------------------------------|
167 subroutine dgesl
(a
,lda
,n
,ipvt
,b
,job
)
168 integer lda
,n
,ipvt
(n
),job
169 double precision a
(lda
,n
),b
(n
)
171 c dgesl solves the double precision system
172 c a * x = b or trans(a) * x = b
173 c using the factors computed by dgeco or dgefa.
177 c a double precision(lda, n)
178 c the output from dgeco or dgefa.
181 c the leading dimension of the array a .
184 c the order of the matrix a .
187 c the pivot vector from dgeco or dgefa.
189 c b double precision(n)
190 c the right hand side vector.
193 c = 0 to solve a*x = b ,
194 c = nonzero to solve trans(a)*x = b where
195 c trans(a) is the transpose.
199 c b the solution vector x .
203 c a division by zero will occur if the input factor contains a
204 c zero on the diagonal. technically this indicates singularity
205 c but it is often caused by improper arguments or improper
206 c setting of lda . it will not occur if the subroutines are
207 c called correctly and if dgeco has set rcond .gt. 0.0
208 c or dgefa has set info .eq. 0 .
210 c to compute inverse(a) * c where c is a matrix
212 c call dgeco(a,lda,n,ipvt,rcond,z)
213 c if (rcond is too small) go to ...
215 c call dgesl(a,lda,n,ipvt,c(1,j),0)
218 c linpack. this version dated 08/14/78 .
219 c cleve moler, university of new mexico, argonne national lab.
221 c subroutines and functions
227 double precision ddot
,t
231 if (job
.ne
. 0) go to 50
233 c job = 0 , solve a * x = b
234 c first solve l*y = b
236 if (nm1
.lt
. 1) go to 30
240 if (l
.eq
. k
) go to 10
244 call daxpy
(n
-k
,t
,a
(k
+1,k
),1,b
(k
+1),1)
254 call daxpy
(k
-1,t
,a
(1,k
),1,b
(1),1)
259 c job = nonzero, solve trans(a) * x = b
260 c first solve trans(u)*y = b
263 t
= ddot
(k
-1,a
(1,k
),1,b
(1),1)
264 b
(k
) = (b
(k
) - t
)/a
(k
,k
)
267 c now solve trans(l)*x = y
269 if (nm1
.lt
. 1) go to 90
272 b
(k
) = b
(k
) + ddot
(n
-k
,a
(k
+1,k
),1,b
(k
+1),1)
274 if (l
.eq
. k
) go to 70
285 subroutine dsifa
(a
,lda
,n
,kpvt
,info
)
286 integer lda
,n
,kpvt
(1),info
287 double precision a
(lda
,1)
289 c dsifa factors a double precision symmetric matrix by elimination
290 c with symmetric pivoting.
292 c to solve a*x = b , follow dsifa by dsisl.
293 c to compute inverse(a)*c , follow dsifa by dsisl.
294 c to compute determinant(a) , follow dsifa by dsidi.
295 c to compute inertia(a) , follow dsifa by dsidi.
296 c to compute inverse(a) , follow dsifa by dsidi.
300 c a double precision(lda,n)
301 c the symmetric matrix to be factored.
302 c only the diagonal and upper triangle are used.
305 c the leading dimension of the array a .
308 c the order of the matrix a .
312 c a a block diagonal matrix and the multipliers which
313 c were used to obtain it.
314 c the factorization can be written a = u*d*trans(u)
315 c where u is a product of permutation and unit
316 c upper triangular matrices , trans(u) is the
317 c transpose of u , and d is block diagonal
318 c with 1 by 1 and 2 by 2 blocks.
321 c an integer vector of pivot indices.
325 c = k if the k-th pivot block is singular. this is
326 c not an error condition for this subroutine,
327 c but it does indicate that dsisl or dsidi may
328 c divide by zero if called.
330 c linpack. this version dated 08/14/78 .
331 c james bunch, univ. calif. san diego, argonne nat. lab.
333 c subroutines and functions
335 c blas daxpy,dswap,idamax
336 c fortran dabs,dmax1,dsqrt
340 double precision ak
,akm1
,bk
,bkm1
,denom
,mulk
,mulkm1
,t
341 double precision absakk
,alpha
,colmax
,rowmax
342 integer imax
,imaxp1
,j
,jj
,jmax
,k
,km1
,km2
,kstep
,idamax
348 c alpha is used in choosing pivot block size.
349 alpha
= (1.0d0
+ dsqrt
(17.0d0
))/8.0d0
353 c main loop on k, which goes from n to 1.
358 c leave the loop if k=0 or k=1.
361 if (k
.eq
. 0) go to 200
362 if (k
.gt
. 1) go to 20
364 if (a
(1,1) .eq
. 0.0d0
) info
= 1
369 c this section of code determines the kind of
370 c elimination to be performed. when it is completed,
371 c kstep will be set to the size of the pivot block, and
372 c swap will be set to .true. if an interchange is
376 absakk
= dabs
(a
(k
,k
))
378 c determine the largest off-diagonal element in
381 imax
= idamax
(k
-1,a
(1,k
),1)
382 colmax
= dabs
(a
(imax
,k
))
383 if (absakk
.lt
. alpha*colmax
) go to 30
389 c determine the largest off-diagonal element in
395 rowmax
= dmax1
(rowmax
,dabs
(a
(imax
,j
)))
397 if (imax
.eq
. 1) go to 50
398 jmax
= idamax
(imax
-1,a
(1,imax
),1)
399 rowmax
= dmax1
(rowmax
,dabs
(a
(jmax
,imax
)))
401 if (dabs
(a
(imax
,imax
)) .lt
. alpha*rowmax
) go to 60
406 if (absakk
.lt
. alpha*colmax*
(colmax
/rowmax
)) go to 70
415 if (dmax1
(absakk
,colmax
) .ne
. 0.0d0
) go to 100
417 c column k is zero. set info and iterate the loop.
423 if (kstep
.eq
. 2) go to 140
427 if (.not
.swap
) go to 120
429 c perform an interchange.
431 call dswap
(imax
,a
(1,imax
),1,a
(1,k
),1)
440 c perform the elimination.
444 mulk
= -a
(j
,k
)/a
(k
,k
)
446 call daxpy
(j
,t
,a
(1,k
),1,a
(1,j
),1)
450 c set the pivot array.
453 if (swap
) kpvt
(k
) = imax
459 if (.not
.swap
) go to 160
461 c perform an interchange.
463 call dswap
(imax
,a
(1,imax
),1,a
(1,k
-1),1)
464 do 150 jj
= imax
, km1
475 c perform the elimination.
478 if (km2
.eq
. 0) go to 180
480 akm1
= a
(k
-1,k
-1)/a
(k
-1,k
)
481 denom
= 1.0d0
- ak*akm1
485 bkm1
= a
(j
,k
-1)/a
(k
-1,k
)
486 mulk
= (akm1*bk
- bkm1
)/denom
487 mulkm1
= (ak*bkm1
- bk
)/denom
489 call daxpy
(j
,t
,a
(1,k
),1,a
(1,j
),1)
491 call daxpy
(j
,t
,a
(1,k
-1),1,a
(1,j
),1)
497 c set the pivot array.
500 if (swap
) kpvt
(k
) = -imax
508 subroutine dsisl
(a
,lda
,n
,kpvt
,b
)
509 integer lda
,n
,kpvt
(1)
510 double precision a
(lda
,1),b
(1)
512 c dsisl solves the double precision symmetric system
514 c using the factors computed by dsifa.
518 c a double precision(lda,n)
519 c the output from dsifa.
522 c the leading dimension of the array a .
525 c the order of the matrix a .
528 c the pivot vector from dsifa.
530 c b double precision(n)
531 c the right hand side vector.
535 c b the solution vector x .
539 c a division by zero may occur if dsico has set rcond .eq. 0.0
540 c or dsifa has set info .ne. 0 .
542 c to compute inverse(a) * c where c is a matrix
544 c call dsifa(a,lda,n,kpvt,info)
545 c if (info .ne. 0) go to ...
547 c call dsisl(a,lda,n,kpvt,c(1,j))
550 c linpack. this version dated 08/14/78 .
551 c james bunch, univ. calif. san diego, argonne nat. lab.
553 c subroutines and functions
558 c internal variables.
560 double precision ak
,akm1
,bk
,bkm1
,ddot
,denom
,temp
563 c loop backward applying the transformations and
567 10 if (k
.eq
. 0) go to 80
568 if (kpvt
(k
) .lt
. 0) go to 40
572 if (k
.eq
. 1) go to 30
574 if (kp
.eq
. k
) go to 20
583 c apply the transformation.
585 call daxpy
(k
-1,b
(k
),a
(1,k
),1,b
(1),1)
597 if (k
.eq
. 2) go to 60
599 if (kp
.eq
. k
- 1) go to 50
608 c apply the transformation.
610 call daxpy
(k
-2,b
(k
),a
(1,k
),1,b
(1),1)
611 call daxpy
(k
-2,b
(k
-1),a
(1,k
-1),1,b
(1),1)
617 akm1
= a
(k
-1,k
-1)/a
(k
-1,k
)
619 bkm1
= b
(k
-1)/a
(k
-1,k
)
620 denom
= ak*akm1
- 1.0d0
621 b
(k
) = (akm1*bk
- bkm1
)/denom
622 b
(k
-1) = (ak*bkm1
- bk
)/denom
628 c loop forward applying the transformations.
631 90 if (k
.gt
. n
) go to 160
632 if (kpvt
(k
) .lt
. 0) go to 120
636 if (k
.eq
. 1) go to 110
638 c apply the transformation.
640 b
(k
) = b
(k
) + ddot
(k
-1,a
(1,k
),1,b
(1),1)
642 if (kp
.eq
. k
) go to 100
657 if (k
.eq
. 1) go to 140
659 c apply the transformation.
661 b
(k
) = b
(k
) + ddot
(k
-1,a
(1,k
),1,b
(1),1)
662 b
(k
+1) = b
(k
+1) + ddot
(k
-1,a
(1,k
+1),1,b
(1),1)
664 if (kp
.eq
. k
) go to 130
679 subroutine zgefa
(a
,lda
,n
,ipvt
,info
)
680 integer lda
,n
,ipvt
(1),info
683 c zgefa factors a complex*16 matrix by gaussian elimination.
685 c zgefa is usually called by zgeco, but it can be called
686 c directly with a saving in time if rcond is not needed.
687 c (time for zgeco) = (1 + 9/n)*(time for zgefa) .
691 c a complex*16(lda, n)
692 c the matrix to be factored.
695 c the leading dimension of the array a .
698 c the order of the matrix a .
702 c a an upper triangular matrix and the multipliers
703 c which were used to obtain it.
704 c the factorization can be written a = l*u where
705 c l is a product of permutation and unit lower
706 c triangular matrices and u is upper triangular.
709 c an integer vector of pivot indices.
713 c = k if u(k,k) .eq. 0.0 . this is not an error
714 c condition for this subroutine, but it does
715 c indicate that zgesl or zgedi will divide by zero
716 c if called. use rcond in zgeco for a reliable
717 c indication of singularity.
719 c linpack. this version dated 08/14/78 .
720 c cleve moler, university of new mexico, argonne national lab.
722 c subroutines and functions
724 c blas zaxpy,zscal,izamax
730 integer izamax
,j
,k
,kp1
,l
,nm1
733 double precision cabs1
734 double precision dreal
,dimag
735 complex*16 zdumr
,zdumi
737 dimag
(zdumi
) = (0.0d0
,-1.0d0
)*zdumi
738 cabs1
(zdum
) = dabs
(dreal
(zdum
)) + dabs
(dimag
(zdum
))
740 c gaussian elimination with partial pivoting
744 if (nm1
.lt
. 1) go to 70
748 c find l = pivot index
750 l
= izamax
(n
-k
+1,a
(k
,k
),1) + k
- 1
753 c zero pivot implies this column already triangularized
755 if (cabs1
(a
(l
,k
)) .eq
. 0.0d0
) go to 40
757 c interchange if necessary
759 if (l
.eq
. k
) go to 10
765 c compute multipliers
767 t
= -(1.0d0
,0.0d0
)/a
(k
,k
)
768 call zscal
(n
-k
,t
,a
(k
+1,k
),1)
770 c row elimination with column indexing
774 if (l
.eq
. k
) go to 20
778 call zaxpy
(n
-k
,t
,a
(k
+1,k
),1,a
(k
+1,j
),1)
787 if (cabs1
(a
(n
,n
)) .eq
. 0.0d0
) info
= n
790 subroutine zgesl
(a
,lda
,n
,ipvt
,b
,job
)
791 integer lda
,n
,ipvt
(1),job
792 complex*16 a
(lda
,1),b
(1)
794 c zgesl solves the complex*16 system
795 c a * x = b or ctrans(a) * x = b
796 c using the factors computed by zgeco or zgefa.
800 c a complex*16(lda, n)
801 c the output from zgeco or zgefa.
804 c the leading dimension of the array a .
807 c the order of the matrix a .
810 c the pivot vector from zgeco or zgefa.
813 c the right hand side vector.
816 c = 0 to solve a*x = b ,
817 c = nonzero to solve ctrans(a)*x = b where
818 c ctrans(a) is the conjugate transpose.
822 c b the solution vector x .
826 c a division by zero will occur if the input factor contains a
827 c zero on the diagonal. technically this indicates singularity
828 c but it is often caused by improper arguments or improper
829 c setting of lda . it will not occur if the subroutines are
830 c called correctly and if zgeco has set rcond .gt. 0.0
831 c or zgefa has set info .eq. 0 .
833 c to compute inverse(a) * c where c is a matrix
835 c call zgeco(a,lda,n,ipvt,rcond,z)
836 c if (rcond is too small) go to ...
838 c call zgesl(a,lda,n,ipvt,c(1,j),0)
841 c linpack. this version dated 08/14/78 .
842 c cleve moler, university of new mexico, argonne national lab.
844 c subroutines and functions
853 double precision dreal
,dimag
854 complex*16 zdumr
,zdumi
856 dimag
(zdumi
) = (0.0d0
,-1.0d0
)*zdumi
859 if (job
.ne
. 0) go to 50
861 c job = 0 , solve a * x = b
862 c first solve l*y = b
864 if (nm1
.lt
. 1) go to 30
868 if (l
.eq
. k
) go to 10
872 call zaxpy
(n
-k
,t
,a
(k
+1,k
),1,b
(k
+1),1)
882 call zaxpy
(k
-1,t
,a
(1,k
),1,b
(1),1)
887 c job = nonzero, solve ctrans(a) * x = b
888 c first solve ctrans(u)*y = b
891 t
= zdotc
(k
-1,a
(1,k
),1,b
(1),1)
892 b
(k
) = (b
(k
) - t
)/dconjg
(a
(k
,k
))
895 c now solve ctrans(l)*x = y
897 if (nm1
.lt
. 1) go to 90
900 b
(k
) = b
(k
) + zdotc
(n
-k
,a
(k
+1,k
),1,b
(k
+1),1)
902 if (l
.eq
. k
) go to 70
912 subroutine zhifa
(a
,lda
,n
,kpvt
,info
)
913 integer lda
,n
,kpvt
(1),info
916 c zhifa factors a complex*16 hermitian matrix by elimination
917 c with symmetric pivoting.
919 c to solve a*x = b , follow zhifa by zhisl.
920 c to compute inverse(a)*c , follow zhifa by zhisl.
921 c to compute determinant(a) , follow zhifa by zhidi.
922 c to compute inertia(a) , follow zhifa by zhidi.
923 c to compute inverse(a) , follow zhifa by zhidi.
927 c a complex*16(lda,n)
928 c the hermitian matrix to be factored.
929 c only the diagonal and upper triangle are used.
932 c the leading dimension of the array a .
935 c the order of the matrix a .
939 c a a block diagonal matrix and the multipliers which
940 c were used to obtain it.
941 c the factorization can be written a = u*d*ctrans(u)
942 c where u is a product of permutation and unit
943 c upper triangular matrices , ctrans(u) is the
944 c conjugate transpose of u , and d is block diagonal
945 c with 1 by 1 and 2 by 2 blocks.
948 c an integer vector of pivot indices.
952 c = k if the k-th pivot block is singular. this is
953 c not an error condition for this subroutine,
954 c but it does indicate that zhisl or zhidi may
955 c divide by zero if called.
957 c linpack. this version dated 08/14/78 .
958 c james bunch, univ. calif. san diego, argonne nat. lab.
960 c subroutines and functions
962 c blas zaxpy,zswap,izamax
963 c fortran dabs,dmax1,dcmplx,dconjg,dsqrt
967 complex*16 ak
,akm1
,bk
,bkm1
,denom
,mulk
,mulkm1
,t
968 double precision absakk
,alpha
,colmax
,rowmax
969 integer imax
,imaxp1
,j
,jj
,jmax
,k
,km1
,km2
,kstep
,izamax
973 double precision cabs1
974 double precision dreal
,dimag
975 complex*16 zdumr
,zdumi
977 dimag
(zdumi
) = (0.0d0
,-1.0d0
)*zdumi
978 cabs1
(zdum
) = dabs
(dreal
(zdum
)) + dabs
(dimag
(zdum
))
982 c alpha is used in choosing pivot block size.
983 alpha
= (1.0d0
+ dsqrt
(17.0d0
))/8.0d0
987 c main loop on k, which goes from n to 1.
992 c leave the loop if k=0 or k=1.
995 if (k
.eq
. 0) go to 200
996 if (k
.gt
. 1) go to 20
998 if (cabs1
(a
(1,1)) .eq
. 0.0d0
) info
= 1
1003 c this section of code determines the kind of
1004 c elimination to be performed. when it is completed,
1005 c kstep will be set to the size of the pivot block, and
1006 c swap will be set to .true. if an interchange is
1010 absakk
= cabs1
(a
(k
,k
))
1012 c determine the largest off-diagonal element in
1015 imax
= izamax
(k
-1,a
(1,k
),1)
1016 colmax
= cabs1
(a
(imax
,k
))
1017 if (absakk
.lt
. alpha*colmax
) go to 30
1023 c determine the largest off-diagonal element in
1029 rowmax
= dmax1
(rowmax
,cabs1
(a
(imax
,j
)))
1031 if (imax
.eq
. 1) go to 50
1032 jmax
= izamax
(imax
-1,a
(1,imax
),1)
1033 rowmax
= dmax1
(rowmax
,cabs1
(a
(jmax
,imax
)))
1035 if (cabs1
(a
(imax
,imax
)) .lt
. alpha*rowmax
) go to 60
1040 if (absakk
.lt
. alpha*colmax*
(colmax
/rowmax
)) go to 70
1046 swap
= imax
.ne
. km1
1049 if (dmax1
(absakk
,colmax
) .ne
. 0.0d0
) go to 100
1051 c column k is zero. set info and iterate the loop.
1057 if (kstep
.eq
. 2) go to 140
1059 c 1 x 1 pivot block.
1061 if (.not
.swap
) go to 120
1063 c perform an interchange.
1065 call zswap
(imax
,a
(1,imax
),1,a
(1,k
),1)
1069 a
(j
,k
) = dconjg
(a
(imax
,j
))
1074 c perform the elimination.
1078 mulk
= -a
(j
,k
)/a
(k
,k
)
1080 call zaxpy
(j
,t
,a
(1,k
),1,a
(1,j
),1)
1081 a
(j
,j
) = dcmplx
(dreal
(a
(j
,j
)),0.0d0
)
1085 c set the pivot array.
1088 if (swap
) kpvt
(k
) = imax
1092 c 2 x 2 pivot block.
1094 if (.not
.swap
) go to 160
1096 c perform an interchange.
1098 call zswap
(imax
,a
(1,imax
),1,a
(1,k
-1),1)
1099 do 150 jj
= imax
, km1
1101 t
= dconjg
(a
(j
,k
-1))
1102 a
(j
,k
-1) = dconjg
(a
(imax
,j
))
1106 a
(k
-1,k
) = a
(imax
,k
)
1110 c perform the elimination.
1113 if (km2
.eq
. 0) go to 180
1114 ak
= a
(k
,k
)/a
(k
-1,k
)
1115 akm1
= a
(k
-1,k
-1)/dconjg
(a
(k
-1,k
))
1116 denom
= 1.0d0
- ak*akm1
1119 bk
= a
(j
,k
)/a
(k
-1,k
)
1120 bkm1
= a
(j
,k
-1)/dconjg
(a
(k
-1,k
))
1121 mulk
= (akm1*bk
- bkm1
)/denom
1122 mulkm1
= (ak*bkm1
- bk
)/denom
1124 call zaxpy
(j
,t
,a
(1,k
),1,a
(1,j
),1)
1126 call zaxpy
(j
,t
,a
(1,k
-1),1,a
(1,j
),1)
1129 a
(j
,j
) = dcmplx
(dreal
(a
(j
,j
)),0.0d0
)
1133 c set the pivot array.
1136 if (swap
) kpvt
(k
) = -imax
1144 subroutine zhisl
(a
,lda
,n
,kpvt
,b
)
1145 integer lda
,n
,kpvt
(1)
1146 complex*16 a
(lda
,1),b
(1)
1148 c zhisl solves the complex*16 hermitian system
1150 c using the factors computed by zhifa.
1154 c a complex*16(lda,n)
1155 c the output from zhifa.
1158 c the leading dimension of the array a .
1161 c the order of the matrix a .
1164 c the pivot vector from zhifa.
1167 c the right hand side vector.
1171 c b the solution vector x .
1175 c a division by zero may occur if zhico has set rcond .eq. 0.0
1176 c or zhifa has set info .ne. 0 .
1178 c to compute inverse(a) * c where c is a matrix
1180 c call zhifa(a,lda,n,kpvt,info)
1181 c if (info .ne. 0) go to ...
1183 c call zhisl(a,lda,n,kpvt,c(1,j))
1186 c linpack. this version dated 08/14/78 .
1187 c james bunch, univ. calif. san diego, argonne nat. lab.
1189 c subroutines and functions
1192 c fortran dconjg,iabs
1194 c internal variables.
1196 complex*16 ak
,akm1
,bk
,bkm1
,zdotc
,denom
,temp
1199 c loop backward applying the transformations and
1203 10 if (k
.eq
. 0) go to 80
1204 if (kpvt
(k
) .lt
. 0) go to 40
1206 c 1 x 1 pivot block.
1208 if (k
.eq
. 1) go to 30
1210 if (kp
.eq
. k
) go to 20
1219 c apply the transformation.
1221 call zaxpy
(k
-1,b
(k
),a
(1,k
),1,b
(1),1)
1231 c 2 x 2 pivot block.
1233 if (k
.eq
. 2) go to 60
1235 if (kp
.eq
. k
- 1) go to 50
1244 c apply the transformation.
1246 call zaxpy
(k
-2,b
(k
),a
(1,k
),1,b
(1),1)
1247 call zaxpy
(k
-2,b
(k
-1),a
(1,k
-1),1,b
(1),1)
1252 ak
= a
(k
,k
)/dconjg
(a
(k
-1,k
))
1253 akm1
= a
(k
-1,k
-1)/a
(k
-1,k
)
1254 bk
= b
(k
)/dconjg
(a
(k
-1,k
))
1255 bkm1
= b
(k
-1)/a
(k
-1,k
)
1256 denom
= ak*akm1
- 1.0d0
1257 b
(k
) = (akm1*bk
- bkm1
)/denom
1258 b
(k
-1) = (ak*bkm1
- bk
)/denom
1264 c loop forward applying the transformations.
1267 90 if (k
.gt
. n
) go to 160
1268 if (kpvt
(k
) .lt
. 0) go to 120
1270 c 1 x 1 pivot block.
1272 if (k
.eq
. 1) go to 110
1274 c apply the transformation.
1276 b
(k
) = b
(k
) + zdotc
(k
-1,a
(1,k
),1,b
(1),1)
1278 if (kp
.eq
. k
) go to 100
1291 c 2 x 2 pivot block.
1293 if (k
.eq
. 1) go to 140
1295 c apply the transformation.
1297 b
(k
) = b
(k
) + zdotc
(k
-1,a
(1,k
),1,b
(1),1)
1298 b
(k
+1) = b
(k
+1) + zdotc
(k
-1,a
(1,k
+1),1,b
(1),1)
1300 if (kp
.eq
. k
) go to 130
1315 subroutine zsifa
(a
,lda
,n
,kpvt
,info
)
1316 integer lda
,n
,kpvt
(1),info
1319 c zsifa factors a complex*16 symmetric matrix by elimination
1320 c with symmetric pivoting.
1322 c to solve a*x = b , follow zsifa by zsisl.
1323 c to compute inverse(a)*c , follow zsifa by zsisl.
1324 c to compute determinant(a) , follow zsifa by zsidi.
1325 c to compute inverse(a) , follow zsifa by zsidi.
1329 c a complex*16(lda,n)
1330 c the symmetric matrix to be factored.
1331 c only the diagonal and upper triangle are used.
1334 c the leading dimension of the array a .
1337 c the order of the matrix a .
1341 c a a block diagonal matrix and the multipliers which
1342 c were used to obtain it.
1343 c the factorization can be written a = u*d*trans(u)
1344 c where u is a product of permutation and unit
1345 c upper triangular matrices , trans(u) is the
1346 c transpose of u , and d is block diagonal
1347 c with 1 by 1 and 2 by 2 blocks.
1350 c an integer vector of pivot indices.
1354 c = k if the k-th pivot block is singular. this is
1355 c not an error condition for this subroutine,
1356 c but it does indicate that zsisl or zsidi may
1357 c divide by zero if called.
1359 c linpack. this version dated 08/14/78 .
1360 c james bunch, univ. calif. san diego, argonne nat. lab.
1362 c subroutines and functions
1364 c blas zaxpy,zswap,izamax
1365 c fortran dabs,dmax1,dsqrt
1367 c internal variables
1369 complex*16 ak
,akm1
,bk
,bkm1
,denom
,mulk
,mulkm1
,t
1370 double precision absakk
,alpha
,colmax
,rowmax
1371 integer imax
,imaxp1
,j
,jj
,jmax
,k
,km1
,km2
,kstep
,izamax
1375 double precision cabs1
1376 double precision dreal
,dimag
1377 complex*16 zdumr
,zdumi
1378 dreal
(zdumr
) = zdumr
1379 dimag
(zdumi
) = (0.0d0
,-1.0d0
)*zdumi
1380 cabs1
(zdum
) = dabs
(dreal
(zdum
)) + dabs
(dimag
(zdum
))
1384 c alpha is used in choosing pivot block size.
1385 alpha
= (1.0d0
+ dsqrt
(17.0d0
))/8.0d0
1389 c main loop on k, which goes from n to 1.
1394 c leave the loop if k=0 or k=1.
1397 if (k
.eq
. 0) go to 200
1398 if (k
.gt
. 1) go to 20
1400 if (cabs1
(a
(1,1)) .eq
. 0.0d0
) info
= 1
1405 c this section of code determines the kind of
1406 c elimination to be performed. when it is completed,
1407 c kstep will be set to the size of the pivot block, and
1408 c swap will be set to .true. if an interchange is
1412 absakk
= cabs1
(a
(k
,k
))
1414 c determine the largest off-diagonal element in
1417 imax
= izamax
(k
-1,a
(1,k
),1)
1418 colmax
= cabs1
(a
(imax
,k
))
1419 if (absakk
.lt
. alpha*colmax
) go to 30
1425 c determine the largest off-diagonal element in
1431 rowmax
= dmax1
(rowmax
,cabs1
(a
(imax
,j
)))
1433 if (imax
.eq
. 1) go to 50
1434 jmax
= izamax
(imax
-1,a
(1,imax
),1)
1435 rowmax
= dmax1
(rowmax
,cabs1
(a
(jmax
,imax
)))
1437 if (cabs1
(a
(imax
,imax
)) .lt
. alpha*rowmax
) go to 60
1442 if (absakk
.lt
. alpha*colmax*
(colmax
/rowmax
)) go to 70
1448 swap
= imax
.ne
. km1
1451 if (dmax1
(absakk
,colmax
) .ne
. 0.0d0
) go to 100
1453 c column k is zero. set info and iterate the loop.
1459 if (kstep
.eq
. 2) go to 140
1461 c 1 x 1 pivot block.
1463 if (.not
.swap
) go to 120
1465 c perform an interchange.
1467 call zswap
(imax
,a
(1,imax
),1,a
(1,k
),1)
1476 c perform the elimination.
1480 mulk
= -a
(j
,k
)/a
(k
,k
)
1482 call zaxpy
(j
,t
,a
(1,k
),1,a
(1,j
),1)
1486 c set the pivot array.
1489 if (swap
) kpvt
(k
) = imax
1493 c 2 x 2 pivot block.
1495 if (.not
.swap
) go to 160
1497 c perform an interchange.
1499 call zswap
(imax
,a
(1,imax
),1,a
(1,k
-1),1)
1500 do 150 jj
= imax
, km1
1503 a
(j
,k
-1) = a
(imax
,j
)
1507 a
(k
-1,k
) = a
(imax
,k
)
1511 c perform the elimination.
1514 if (km2
.eq
. 0) go to 180
1515 ak
= a
(k
,k
)/a
(k
-1,k
)
1516 akm1
= a
(k
-1,k
-1)/a
(k
-1,k
)
1517 denom
= 1.0d0
- ak*akm1
1520 bk
= a
(j
,k
)/a
(k
-1,k
)
1521 bkm1
= a
(j
,k
-1)/a
(k
-1,k
)
1522 mulk
= (akm1*bk
- bkm1
)/denom
1523 mulkm1
= (ak*bkm1
- bk
)/denom
1525 call zaxpy
(j
,t
,a
(1,k
),1,a
(1,j
),1)
1527 call zaxpy
(j
,t
,a
(1,k
-1),1,a
(1,j
),1)
1533 c set the pivot array.
1536 if (swap
) kpvt
(k
) = -imax
1544 subroutine zsisl
(a
,lda
,n
,kpvt
,b
)
1545 integer lda
,n
,kpvt
(1)
1546 complex*16 a
(lda
,1),b
(1)
1548 c zsisl solves the complex*16 symmetric system
1550 c using the factors computed by zsifa.
1554 c a complex*16(lda,n)
1555 c the output from zsifa.
1558 c the leading dimension of the array a .
1561 c the order of the matrix a .
1564 c the pivot vector from zsifa.
1567 c the right hand side vector.
1571 c b the solution vector x .
1575 c a division by zero may occur if zsico has set rcond .eq. 0.0
1576 c or zsifa has set info .ne. 0 .
1578 c to compute inverse(a) * c where c is a matrix
1580 c call zsifa(a,lda,n,kpvt,info)
1581 c if (info .ne. 0) go to ...
1583 c call zsisl(a,lda,n,kpvt,c(1,j))
1586 c linpack. this version dated 08/14/78 .
1587 c james bunch, univ. calif. san diego, argonne nat. lab.
1589 c subroutines and functions
1594 c internal variables.
1596 complex*16 ak
,akm1
,bk
,bkm1
,zdotu
,denom
,temp
1599 c loop backward applying the transformations and
1603 10 if (k
.eq
. 0) go to 80
1604 if (kpvt
(k
) .lt
. 0) go to 40
1606 c 1 x 1 pivot block.
1608 if (k
.eq
. 1) go to 30
1610 if (kp
.eq
. k
) go to 20
1619 c apply the transformation.
1621 call zaxpy
(k
-1,b
(k
),a
(1,k
),1,b
(1),1)
1631 c 2 x 2 pivot block.
1633 if (k
.eq
. 2) go to 60
1635 if (kp
.eq
. k
- 1) go to 50
1644 c apply the transformation.
1646 call zaxpy
(k
-2,b
(k
),a
(1,k
),1,b
(1),1)
1647 call zaxpy
(k
-2,b
(k
-1),a
(1,k
-1),1,b
(1),1)
1652 ak
= a
(k
,k
)/a
(k
-1,k
)
1653 akm1
= a
(k
-1,k
-1)/a
(k
-1,k
)
1655 bkm1
= b
(k
-1)/a
(k
-1,k
)
1656 denom
= ak*akm1
- 1.0d0
1657 b
(k
) = (akm1*bk
- bkm1
)/denom
1658 b
(k
-1) = (ak*bkm1
- bk
)/denom
1664 c loop forward applying the transformations.
1667 90 if (k
.gt
. n
) go to 160
1668 if (kpvt
(k
) .lt
. 0) go to 120
1670 c 1 x 1 pivot block.
1672 if (k
.eq
. 1) go to 110
1674 c apply the transformation.
1676 b
(k
) = b
(k
) + zdotu
(k
-1,a
(1,k
),1,b
(1),1)
1678 if (kp
.eq
. k
) go to 100
1691 c 2 x 2 pivot block.
1693 if (k
.eq
. 1) go to 140
1695 c apply the transformation.
1697 b
(k
) = b
(k
) + zdotu
(k
-1,a
(1,k
),1,b
(1),1)
1698 b
(k
+1) = b
(k
+1) + zdotu
(k
-1,a
(1,k
+1),1,b
(1),1)
1700 if (kp
.eq
. k
) go to 130