MODULE LinAlg USE IntrType ! ! This module contains some very old implementations of LINPACK ! LU factorization/solution subroutines without reference to BLAS ! Also contains a simple tridiagonal equation solver ! CONTAINS SUBROUTINE dgefa(a,lda,n,ipvt,info) IMPLICIT NONE ! ! Declaration Generated by genImpDecs.pl 5/98 INTEGER(sik) i,info,ip1,j,k,lda,m,n,nm1 ! ! Declaration Generated by genImpDecs.pl 5/98 REAL(sdk) c,rp,t ! ! ! because its name is same as routine in Cray system library. ! ! ! dgefa factors a real matrix by gaussian elimination. ! ! dgefa is usually called by dgeco, but it can be called ! directly with a saving in time if rcond is not needed. ! (time for dgeco) = (1 + 9/n)*(time for dgefa) . ! ! on entry ! ! a real(lda, n) ! the matrix to be factored. ! ! lda integer ! the leading dimension of the array a . ! ! n integer ! the order of the matrix a . ! ! on return ! ! a an upper triangular matrix and the multipliers ! which were used to obtain it. ! the factorization can be written a = l*u where ! l is a product of permutation and unit lower ! triangular matrices and u is upper triangular. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 0 normal value. ! = k if u(k,k) .eq. 0.0 . this is not an error ! condition for this subroutine, but it does ! indicate that dgesl or dgedit will divide by zero ! if called. use rcond in dgeco for a reliable ! indication of singularity. ! ! this version is written in cdc compass for the ftn compiler ! communication convention. results may vary slightly from those ! obtained with a fortran version, and such a version is available ! upon request. this program simulates the fortran sequence.. ! ! subroutine dgefa (a, lda, n, ipvt, info) ! INTEGER(sik), INTENT(OUT) :: ipvt(:) ! INTEGER(sik) ipvt(n) ! REAL(sdk), INTENT(INOUT) :: a(:,:) ! REAL(sdk) a(lda,n) ! info=0 IF (n.GT.1) THEN nm1=n-1 DO i=1,nm1 m=i ip1=i+1 DO k=ip1,n IF (abs(a(k,i)).GT.abs(a(m,i))) m=k ENDDO ipvt(i)=m IF (a(m,i).EQ.0.d0) THEN info=i ELSE IF (m.NE.i) THEN DO k=i,n t=a(i,k) a(i,k)=a(m,k) a(m,k)=t ENDDO ENDIF rp=-1.d0/a(i,i) DO k=ip1,n a(k,i)=a(k,i)*rp ENDDO DO j=ip1,n c=a(i,j) IF (c.NE.0.d0) THEN DO k=ip1,n a(k,j)=a(k,j)+a(k,i)*c ENDDO ENDIF ENDDO ENDIF ENDDO ipvt(n)=n IF (a(n,n).EQ.0.d0) info=n ELSE ipvt(1)=1 IF (a(1,1).EQ.0.d0) info=1 ENDIF RETURN END SUBROUTINE dgefa SUBROUTINE dgesl(a,lda,n,ipvt,b,job) IMPLICIT NONE ! ! Declaration Generated by genImpDecs.pl 5/98 INTEGER(sik) i,ip1,j,jm1,job,k,kb,km1,kp1,l,lda,m,n,nm1 ! ! Declaration Generated by genImpDecs.pl 5/98 REAL(sdk) c,t ! ! ! because its name is same as routine in Cray system library. ! ! ! dgesl solves the real system ! a * x = b or trans(a) * x = b ! using the factors computed by dgeco or dgefa. ! ! on entry ! ! a real(lda, n) ! the output from dgeco or dgefa. ! ! lda integer ! the leading dimension of the array a . ! ! n integer ! the order of the matrix a . ! ! ipvt integer(n) ! the pivot vector from dgeco or dgefa. ! ! b real(n) ! the right hand side vector. ! ! job integer ! = 0 to solve a*x = b , ! = nonzero to solve trans(a)*x = b where ! trans(a) is the transpose. ! ! on return ! ! b the solution vector x . ! ! error condition ! ! a division by zero will occur if the input factor contains a ! zero on the diagonal. technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of lda . it will not occur if the subroutines are ! called correctly and if dgeco has set rcond .gt. 0.0 ! or dgefa has set info .eq. 0 . ! ! to compute inverse(a) * c where c is a matrix ! with p columns ! call dgeco(a,lda,n,ipvt,rcond,z) ! if (rcond is too small) go to ... ! do 10 j = 1, p ! call dgesl(a,lda,n,ipvt,c(1,j),0) ! 10 continue ! ! this version is written in cdc assembly language for the ftn ! compiler communication convention. results may vary slightly ! from those obtained with a fortran version, and such a version ! is available upon request. this program simulates the ! following fortran sequence.. ! ! subroutine dgesl (a, lda, n, ipvt, b, job) ! INTEGER(sik), INTENT(IN) :: ipvt(:) ! INTEGER(sik) ipvt(n) ! REAL(sdk), INTENT(INOUT) :: a(:,:),b(:) ! REAL(sdk) a(lda,n),b(n) ! nm1=n-1 IF (n.LE.1) THEN b(1)=b(1)/a(1,1) ELSEIF (job.EQ.0) THEN DO i=1,nm1 m=ipvt(i) c=b(m) IF (m.NE.i) THEN b(m)=b(i) b(i)=c ENDIF IF (c.NE.0.d0) THEN ip1=i+1 DO k=ip1,n b(k)=b(k)+a(k,i)*c ENDDO ENDIF ENDDO DO i=1,nm1 j=n-i+1 b(j)=b(j)/a(j,j) c=b(j) IF (c.NE.0.d0) THEN jm1=j-1 DO k=1,jm1 b(k)=b(k)-a(k,j)*c ENDDO ENDIF ENDDO b(1)=b(1)/a(1,1) ELSE DO k=1,n t=0.d0 IF (k.NE.1) THEN km1=k-1 DO i=1,km1 t=t+a(i,k)*b(i) ENDDO ENDIF b(k)=(b(k)-t)/a(k,k) ENDDO DO kb=1,nm1 k=n-kb t=0.d0 kp1=k+1 DO i=kp1,n t=t+a(i,k)*b(i) ENDDO b(k)=b(k)+t l=ipvt(k) IF (l.NE.k) THEN t=b(l) b(l)=b(k) b(k)=t ENDIF ENDDO ENDIF RETURN END SUBROUTINE dgesl SUBROUTINE TriSolve(a,n,x) ! ! Solve a tridiagonal system with coefficients stored in "a", and right ! hand side stored in b, ! ! ! Programmed by John Mahaffy PSU/ARL 12/97 ! ! Input ! a - derived type array containing matrix coefficients ! n - number of variables ! x - right hand side of equation ! Output ! ! x - solution ! ! IMPLICIT NONE INTEGER(sik) :: n, i REAL(sdk), DIMENSION (:) :: x REAL(sdk) :: a3, a1, rdiag, xLast REAL(sdk) :: a(:,:) ! ! Initialize Elimination Variables xLast=0.0_sdk ! ! Forward Elimination DO i=1,n rdiag=1.0_sdk/(a(2,i)-a(1,i)*a3) a3=a(3,i)*rdiag x(i)=(x(i)-a(1,i)*xLast )*rdiag a(3,i) = a3 xLast = x(i) ENDDO ! ! back substitution on the tridiagonal systems ! DO i =n-1,1,-1 x(i)=x(i)-a(3,i)*x(i+1) ENDDO END SUBROUTINE TriSolve REAL(sdk) FUNCTION ddot(n,dx,incx,dy,incy) IMPLICIT NONE ! ! ! forms the dot product of two vectors. ! uses unrolled loops for increments equal to one. ! jack dongarra, linpack, 3/11/78. ! REAL(sdk) dtemp REAL(sdk), INTENT(INOUT) :: dx(:),dy(:) INTEGER(sik) i,incx,incy,ix,iy,m,mp1,n ! ddot=0.0d0 dtemp=0.0d0 IF (n.GT.0) THEN IF (incx.EQ.1.AND.incy.EQ.1) THEN ! ! code for both increments equal to 1 ! ! ! clean-up loop ! m=mod(n,5) IF (m.NE.0) THEN DO i=1,m dtemp=dtemp+dx(i)*dy(i) ENDDO IF (n.LT.5) GOTO 60 ENDIF mp1=m+1 DO i=mp1,n,5 dtemp=dtemp+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2)+dx & & (i+3)*dy(i+3)+dx(i+4)*dy(i+4) ENDDO 60 ddot=dtemp ELSE ! ! code for unequal increments or equal increments ! not equal to 1 ! ix=1 iy=1 IF (incx.LT.0) ix=(-n+1)*incx+1 IF (incy.LT.0) iy=(-n+1)*incy+1 DO i=1,n dtemp=dtemp+dx(ix)*dy(iy) ix=ix+incx iy=iy+incy ENDDO ddot=dtemp ENDIF ENDIF RETURN END FUNCTION ddot SUBROUTINE dscal(n,da,dx,incx) IMPLICIT NONE ! ! ! scales a vector by a constant. ! uses unrolled loops for increment equal to one. ! jack dongarra, linpack, 3/11/78. ! REAL(sdk) da,dx(*) INTEGER(sik) i,incx,m,mp1,n,nincx ! IF (n.GT.0) THEN IF (incx.EQ.1) THEN ! ! code for increment equal to 1 ! ! ! clean-up loop ! m=mod(n,5) IF (m.NE.0) THEN DO i=1,m dx(i)=da*dx(i) ENDDO IF (n.LT.5) GOTO 51 ENDIF mp1=m+1 DO i=mp1,n,5 dx(i)=da*dx(i) dx(i+1)=da*dx(i+1) dx(i+2)=da*dx(i+2) dx(i+3)=da*dx(i+3) dx(i+4)=da*dx(i+4) ENDDO ELSE ! ! code for increment not equal to 1 ! nincx=n*incx DO i=1,nincx,incx dx(i)=da*dx(i) ENDDO ENDIF ENDIF 51 RETURN END SUBROUTINE dscal SUBROUTINE daxpy(n,da,dx,incx,dy,incy) IMPLICIT NONE ! ! ! constant times a vector plus a vector. ! uses unrolled loops for increments equal to one. ! jack dongarra, linpack, 3/11/78. ! REAL(sdk), INTENT(INOUT) :: dx(:),dy(:),da ! REAL(sdk) dx(*),dy(*),da INTEGER(sik) i,incx,incy,ix,iy,m,mp1,n ! IF (n.GT.0) THEN IF (da.NE.0.0d0) THEN IF (incx.EQ.1.AND.incy.EQ.1) THEN ! ! code for both increments equal to 1 ! ! ! clean-up loop ! m=mod(n,4) IF (m.NE.0) THEN DO i=1,m dy(i)=dy(i)+da*dx(i) ENDDO IF (n.LT.4) GOTO 51 ENDIF mp1=m+1 DO i=mp1,n,4 dy(i)=dy(i)+da*dx(i) dy(i+1)=dy(i+1)+da*dx(i+1) dy(i+2)=dy(i+2)+da*dx(i+2) dy(i+3)=dy(i+3)+da*dx(i+3) ENDDO ELSE ! ! code for unequal increments or equal increments ! not equal to 1 ! ix=1 iy=1 IF (incx.LT.0) ix=(-n+1)*incx+1 IF (incy.LT.0) iy=(-n+1)*incy+1 DO i=1,n dy(iy)=dy(iy)+da*dx(ix) ix=ix+incx iy=iy+incy ENDDO ENDIF ENDIF ENDIF 51 RETURN END SUBROUTINE daxpy SUBROUTINE dgbfa(abd,lda,n,ml,mu,ipvt,info) IMPLICIT NONE ! !***begin prologue dgbfa !***date written 780814 (yymmdd) !***revision date 861211 (yymmdd) !***category no. d2a2 !***keywords library=slatec(linpack), ! type=double precision(sgbfa-s dgbfa-d cgbfa-c),banded, ! linear algebra,matrix,matrix factorization !***author moler, c. b., (u. of new mexico) !***purpose factors a double precision band matrix by elimination. !***description ! ! dgbfa factors a double precision band matrix by elimination. ! ! dgbfa is usually called by dgbco, but it can be called ! directly with a saving in time if rcond is not needed. ! ! on entry ! ! abd double precision(lda, n) ! contains the matrix in band storage. the columns ! of the matrix are stored in the columns of abd and ! the diagonals of the matrix are stored in rows ! ml+1 through 2*ml+mu+1 of abd . ! see the comments below for details. ! ! lda integer ! the leading dimension of the array abd . ! lda must be .ge. 2*ml + mu + 1 . ! ! n integer ! the order of the original matrix. ! ! ml integer ! number of diagonals below the main diagonal. ! 0 .le. ml .lt. n . ! ! mu integer ! number of diagonals above the main diagonal. ! 0 .le. mu .lt. n . ! more efficient if ml .le. mu . ! on return ! ! abd an upper triangular matrix in band storage and ! the multipliers which were used to obtain it. ! the factorization can be written a = l*u where ! l is a product of permutation and unit lower ! triangular matrices and u is upper triangular. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 0 normal value. ! = k if u(k,k) .eq. 0.0 . this is not an error ! condition for this subroutine, but it does ! indicate that dgbsl will divide by zero if ! called. use rcond in dgbco for a reliable ! indication of singularity. ! ! band storage ! ! if a is a band matrix, the following program segment ! will set up the input. ! ! ml = (band width below the diagonal) ! mu = (band width above the diagonal) ! m = ml + mu + 1 ! do 20 j = 1, n ! i1 = max0(1, j-mu) ! i2 = min0(n, j+ml) ! do 10 i = i1, i2 ! k = i - j + m ! abd(k,j) = a(i,j) ! 10 continue ! 20 continue ! ! this uses rows ml+1 through 2*ml+mu+1 of abd . ! in addition, the first ml rows in abd are used for ! elements generated during the triangularization. ! the total number of rows needed in abd is 2*ml+mu+1 . ! the ml+mu by ml+mu upper left triangle and the ! ml by ml lower right triangle are not referenced. ! ! linpack. this version dated 08/14/78 . ! cleve moler, university of new mexico, argonne national lab. ! ! subroutines and functions ! ! blas daxpy,dscal,idamax ! fortran max0,min0 !***references dongarra j.j., bunch j.r., moler c.b., stewart g.w., ! *linpack users guide*, siam, 1979. !***routines called daxpy,dscal,idamax !***end prologue dgbfa INTEGER(sik) lda,n,ml,mu,ipvt(*),info REAL(sdk) abd(lda,n) ! REAL(sdk) t INTEGER(sik) i,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1 ! !***first executable statement dgbfa m=ml+mu+1 info=0 ! ! zero initial fill-in columns ! j0=mu+2 j1= min(n,m)-1 IF (j1.GE.j0) THEN DO jz=j0,j1 i0=m+1-jz DO i=i0,ml abd(i,jz)=0.0d0 ENDDO ENDDO ENDIF jz=j1 ju=0 ! ! gaussian elimination with partial pivoting ! nm1=n-1 IF (nm1.GE.1) THEN DO k=1,nm1 kp1=k+1 ! ! zero next fill-in column ! jz=jz+1 IF (jz.LE.n) THEN IF (ml.GE.1) THEN DO i=1,ml abd(i,jz)=0.0d0 ENDDO ENDIF ENDIF ! ! find l = pivot index ! lm= min(ml,n-k) l=idamax(lm+1,abd(m,k),1)+m-1 ipvt(k)=l+k-m ! ! zero pivot implies this column already triangularized ! IF (abd(l,k).EQ.0.0d0) THEN info=k ELSE ! ! interchange if necessary ! IF (l.NE.m) THEN t=abd(l,k) abd(l,k)=abd(m,k) abd(m,k)=t ENDIF ! ! compute multipliers ! t=-1.0d0/abd(m,k) CALL dscal(lm,t,abd(m+1,k),1) ! ! row elimination with column indexing ! ju= min( max(ju,mu+ipvt(k)),n) mm=m IF (ju.GE.kp1) THEN DO j=kp1,ju l=l-1 mm=mm-1 t=abd(l,j) IF (l.NE.mm) THEN abd(l,j)=abd(mm,j) abd(mm,j)=t ENDIF CALL daxpy(lm,t,abd(m+1:m+lm,k),1,abd(mm+1:mm+lm,j),1) ENDDO ENDIF ENDIF ENDDO ENDIF ipvt(n)=n IF (abd(m,n).EQ.0.0d0) info=n RETURN END SUBROUTINE dgbfa SUBROUTINE dgbsl(abd,lda,n,ml,mu,ipvt,b,job) IMPLICIT NONE ! !***begin prologue dgbsl !***date written 780814 (yymmdd) !***revision date 861211 (yymmdd) !***category no. d2a2 !***keywords library=slatec(linpack), ! type=double precision(sgbsl-s dgbsl-d cgbsl-c),banded, ! linear algebra,matrix,solve !***author moler, c. b., (u. of new mexico) !***purpose solves the double precision band system a*x=b or ! trans(a)*x=b using the factors computed by dgbco or dgbfa. !***description ! ! dgbsl solves the double precision band system ! a * x = b or trans(a) * x = b ! using the factors computed by dgbco or dgbfa. ! ! on entry ! ! abd double precision(lda, n) ! the output from dgbco or dgbfa. ! ! lda integer ! the leading dimension of the array abd . ! ! n integer ! the order of the original matrix. ! ! ml integer ! number of diagonals below the main diagonal. ! ! mu integer ! number of diagonals above the main diagonal. ! ! ipvt integer(n) ! the pivot vector from dgbco or dgbfa. ! ! b double precision(n) ! the right hand side vector. ! ! job integer ! = 0 to solve a*x = b , ! = nonzero to solve trans(a)*x = b , where ! trans(a) is the transpose. ! ! on return ! ! b the solution vector x . ! ! error condition ! ! a division by zero will occur if the input factor contains a ! zero on the diagonal. technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of lda . it will not occur if the subroutines are ! called correctly and if dgbco has set rcond .gt. 0.0 ! or dgbfa has set info .eq. 0 . ! ! to compute inverse(a) * c where c is a matrix ! with p columns ! call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z) ! if (rcond is too small) go to ... ! do 10 j = 1, p ! call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0) ! 10 continue ! ! linpack. this version dated 08/14/78 . ! cleve moler, university of new mexico, argonne national lab. ! ! subroutines and functions ! ! blas daxpy,ddot ! fortran min0 !***references dongarra j.j., bunch j.r., moler c.b., stewart g.w., ! *linpack users guide*, siam, 1979. !***routines called daxpy,ddot !***end prologue dgbsl INTEGER(sik) lda,n,ml,mu,job INTEGER(sik), DIMENSION(:) :: ipvt REAL(sdk), DIMENSION(:) :: b REAL(sdk), DIMENSION(:,:) :: abd ! REAL(sdk) t INTEGER(sik) k,kb,l,la,lb,lm,m,nm1 !***first executable statement dgbsl m=mu+ml+1 nm1=n-1 IF (job.EQ.0) THEN ! ! job = 0 , solve a * x = b ! first solve l*y = b ! IF (ml.NE.0) THEN IF (nm1.GE.1) THEN DO k=1,nm1 lm= min(ml,n-k) l=ipvt(k) t=b(l) IF (l.NE.k) THEN b(l)=b(k) b(k)=t ENDIF CALL daxpy(lm,t,abd(m+1:,k),1,b(k+1:),1) ENDDO ENDIF ENDIF ! ! now solve u*x = y ! DO kb=1,n k=n+1-kb b(k)=b(k)/abd(m,k) lm= min(k,m)-1 la=m-lm lb=k-lm t=-b(k) CALL daxpy(lm,t,abd(la:,k),1,b(lb:),1) ENDDO ELSE ! ! job = nonzero, solve trans(a) * x = b ! first solve trans(u)*y = b ! DO k=1,n lm= min(k,m)-1 la=m-lm lb=k-lm t=ddot(lm,abd(la:,k),1,b(lb:),1) b(k)=(b(k)-t)/abd(m,k) ENDDO ! ! now solve trans(l)*x = y ! IF (ml.NE.0) THEN IF (nm1.GE.1) THEN DO kb=1,nm1 k=n-kb lm= min(ml,n-k) b(k)=b(k)+ddot(lm,abd(m+1:,k),1,b(k+1:),1) l=ipvt(k) IF (l.NE.k) THEN t=b(l) b(l)=b(k) b(k)=t ENDIF ENDDO ENDIF ENDIF ENDIF RETURN END SUBROUTINE dgbsl INTEGER(sik) FUNCTION idamax(n,dx,incx) IMPLICIT NONE ! ! ! finds the index of element having max. absolute value. ! jack dongarra, linpack, 3/11/78. ! REAL(sdk) dx(*),dmax INTEGER(sik) i,incx,ix,n ! idamax=0 IF (n.GE.1) THEN idamax=1 IF (n.NE.1) THEN IF (incx.EQ.1) THEN ! ! code for increment equal to 1 ! dmax=dabs(dx(1)) DO i=2,n IF (dabs(dx(i)).GT.dmax) THEN idamax=i dmax=dabs(dx(i)) ENDIF ENDDO ELSE ! ! code for increment not equal to 1 ! ix=1 dmax=dabs(dx(1)) ix=ix+incx DO i=2,n IF (dabs(dx(ix)).GT.dmax) THEN idamax=i dmax=dabs(dx(ix)) ENDIF ix=ix+incx ENDDO ENDIF ENDIF ENDIF RETURN END FUNCTION idamax ! END MODULE LinAlg