function clock ( ) c*********************************************************************72 c cc clock() returns a reading from the real time clock. c implicit none double precision clock real*4 etime, tm(2) clock = etime( tm ) *--- IBM & CRAY -- * double precision function clock() * real*8 timef * clock = 1000.0d0 * timef() * end *--- others ?? *--- if in trouble, use this to get out of the hook! * double precision function clock() * clock = 0.0d0 * end function dasum ( n, dx, incx ) c*********************************************************************72 c cc DASUM sums the absolute value of the entries of a vector. c c takes the sum of the absolute values. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c implicit none double precision dasum double precision dx(1),dtemp integer i,incx,m,mp1,n,nincx c dasum = 0.0d0 dtemp = 0.0d0 if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dtemp = dtemp + dabs(dx(i)) 10 continue dasum = dtemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end subroutine daxpy ( n, da, dx, incx, dy, incy ) c*********************************************************************72 c cc DAXPY adds a multiple of one vector to another. c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c implicit none double precision dx(1),dy(1),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 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) 50 continue return end function dcabs1 ( z ) c*********************************************************************72 c cc DCABS1 computes the L1 norm of a complex number. c implicit none double precision dcabs1 complex*16 z,zz double precision t(2) equivalence (zz,t(1)) zz = z dcabs1 = dabs(t(1)) + dabs(t(2)) return end subroutine dcmpac ( n, nx, ix, ixx, xx, iwsp ) c*********************************************************************72 c cc DCMPAC compacts the array IX and sorts IXX and XX. c c This is a gateway routine for DGCNVR) ... c implicit none integer n, nx, ix(nx), ixx(nx), iwsp(n) double precision xx(nx) integer k * *--- sort ix and carry ixx and xx along ... * call idsrt2( nx, ix, ixx, xx ) * *--- adjust pointers ... * do k = 1,n iwsp(k) = 0 end do do k = 1,nx iwsp(ix(k)) = iwsp(ix(k)) + 1 end do ix(n+1) = nx + 1 do k = n,1,-1 ix(k) = ix(k+1)-iwsp(k) end do * *--- sort ixx in increasing order and carry xx along ... * do k = 1,n call idsrt1( iwsp(k), ixx(ix(k)), xx(ix(k)) ) end do END subroutine dcopy ( n, dx, incx, dy, incy ) c*********************************************************************72 c cc DCOPY copies a vector. c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c implicit none double precision dx(1),dy(1) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 dy(i) = dx(i) dy(i + 1) = dx(i + 1) dy(i + 2) = dx(i + 2) dy(i + 3) = dx(i + 3) dy(i + 4) = dx(i + 4) dy(i + 5) = dx(i + 5) dy(i + 6) = dx(i + 6) 50 continue return end function ddot(n,dx,incx,dy,incy) c*********************************************************************72 c cc DDOT computes the dot product of two vectors. c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c implicit none double precision ddot double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 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) 50 continue 60 ddot = dtemp return end subroutine dgccsv ( x, y ) c*********************************************************************72 c cc DGCCSV computes y = A * x for A in compressed column storage (CCS). c * *--- Computes y = A*x. A is passed via a fortran `common statement'. *--- A is assumed to be under the Compress Column Storage (CCS) format. * implicit none double precision x(*), y(*) integer n, nz, nzmax parameter( nzmax = 600000 ) integer ia(nzmax), ja(nzmax) double precision a(nzmax) common /RMAT/ a, ia, ja, nz, n integer i, j do i = 1,n y(i) = 0.0d0 end do do j = 1,n do i = ja(j),ja(j+1)-1 y(ia(i)) = y(ia(i)) + a(i)*x(j) end do end do end subroutine dgchbv ( m, t, h,ldh, y, wsp, iwsp, iflag ) c*********************************************************************72 c cc DGCHBV computes y=exp(t*H)*y, general matrix H, using rational Chebyshev. c *-----Purpose----------------------------------------------------------| * *--- DGCHBV computes y = exp(t*H)*y using the partial fraction * expansion of the uniform rational Chebyshev approximation * to exp(-x) of type (14,14). H is a General matrix. * About 14-digit accuracy is expected if the matrix H is negative * definite. The algorithm may behave poorly otherwise. * *-----Arguments--------------------------------------------------------| * * m : (input) order of the matrix H * * t : (input) time-scaling factor (can be < 0). * * H(ldh,m): (input) argument matrix. * * y(m) : (input/output) on input the operand vector, * on output the resulting vector exp(t*H)*y. * * iwsp(m) : (workspace) * * wsp : (workspace). Observe that a double precision vector of * length 2*m*(m+2) can be used as well when calling this * routine (thus avoiding an idle complex array elsewhere) * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * EXPOKIT: Software Package for Computing Matrix Exponentials. * ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 *----------------------------------------------------------------------| * implicit none integer m, ldh, iflag, iwsp(m) double precision t, H(ldh,m), y(m) complex*16 wsp(m*(m+2)) integer ndeg, i, j, ip, ih, iy, iz parameter ( ndeg=7 ) double precision alpha0 complex*16 alpha(ndeg), theta(ndeg) intrinsic DBLE *--- Pointers ... ih = 1 iy = ih + m*m iz = iy + m *--- Coefficients and poles of the partial fraction expansion ... alpha0 = 0.183216998528140087D-11 alpha(1)=( 0.557503973136501826D+02,-0.204295038779771857D+03) alpha(2)=(-0.938666838877006739D+02, 0.912874896775456363D+02) alpha(3)=( 0.469965415550370835D+02,-0.116167609985818103D+02) alpha(4)=(-0.961424200626061065D+01,-0.264195613880262669D+01) alpha(5)=( 0.752722063978321642D+00, 0.670367365566377770D+00) alpha(6)=(-0.188781253158648576D-01,-0.343696176445802414D-01) alpha(7)=( 0.143086431411801849D-03, 0.287221133228814096D-03) theta(1)=(-0.562314417475317895D+01, 0.119406921611247440D+01) theta(2)=(-0.508934679728216110D+01, 0.358882439228376881D+01) theta(3)=(-0.399337136365302569D+01, 0.600483209099604664D+01) theta(4)=(-0.226978543095856366D+01, 0.846173881758693369D+01) theta(5)=( 0.208756929753827868D+00, 0.109912615662209418D+02) theta(6)=( 0.370327340957595652D+01, 0.136563731924991884D+02) theta(7)=( 0.889777151877331107D+01, 0.166309842834712071D+02) * *--- Accumulation of the contribution of each pole ... * do j = 1,m wsp(iz+j-1) = y(j) y(j) = y(j)*alpha0 end do do ip = 1,ndeg *--- Solve each fraction using Gaussian elimination with pivoting... do j = 1,m do i = 1,m wsp(ih+(j-1)*m+i-1) = -t*H(i,j) end do wsp(ih+(j-1)*m+j-1) = wsp(ih+(j-1)*m+j-1)-theta(ip) wsp(iy+j-1) = wsp(iz+j-1) end do call ZGESV( M, 1, WSP(iH),M, IWSP, WSP(iY),M, IFLAG ) if ( IFLAG.ne.0 ) stop 'Error in DGCHBV' c *--- Accumulate the partial result in y. c do j = 1,m y(j) = y(j) + DBLE( alpha(ip)*wsp(iy+j-1) ) end do end do END subroutine dgcnvr ( from, to, diag, nrow, ncol, nz, ia, ja, a, & iwsp ) c*********************************************************************72 c cc DGCNVR converts between sparse matrix storage formats. c *-----Purpose----------------------------------------------------------| * *--- DGCNVR converts a sparse storage format into another sparse * storage format. The matrix can be rectangular. * *-----Arguments--------------------------------------------------------| * * from : (input, character*3) the storage format holding the * matrix on input. Accepted values are: * 'COO' : COOrdinates * 'CRS' : Compressed Row Storage * 'CCS' : Compressed Column Storage. * * to : (input, character*3) the storage format holding the * matrix on output. Same accepted values as above. * * diag : (input, character*1) specifies whether or not the * entire diagonal should be lodged, i.e., including * null diagonal entries. (This may be required by * some applications that make use of the locations * of diagonal elements.) * diag = 'N', no specific attention is given to diagonal * elements. * diag = 'D', the entire diagonal is lodged, including * null elements on the diagonal. * if from=to & diag='D', null diagonal entries are * explicitly inserted in the actual matrix. * * nrow : (input) number of rows in the matrix. * * ncol : (input) number of columns in the matrix. * * nz : (input/output) number of non-zeros elements. * If diag='D' and null diagonal entries are inserted, then * nz is updated on exit and contains the effective number * of entries stored. In what follows, nz' (read nz prime) * denotes the updated value of nz, nz <= nz' <= nz+n. * If diag='N' then nz'=nz. * * ia(*) : (input/output) of declared length .ge. nz'. * On input, * if from='CRS', ia(1:nrow+1) contains pointers for the * beginning of each row. * if from='COO', or 'CCS', ia(1:nz) contains row indices. * On output, * if to='CRS', ia(1:nrow+1) contains pointers for the * beginning of each row. * if to='COO' or 'CCS', ia(1:nz') contains row indices * in increasing order in each column. * * ja(*) : (input/output) of declared length .ge. nz'. * On input, * if from='CRS', ja(1:ncol+1) contains pointers for the * beginning of each column. * if from='COO' or 'CCS', ja(1:nz) contains col. indices. * On output, * if to='CRS', ja(1:ncol+1) contains pointers for the * beginning of each column. * if to='COO' or 'CCS', ja(1:nz') contains column indices * in increasing order in each row. * * a(*) : (input/output) On input, a(1:nz) fits the input format * and on output, a(1:nz') fits the output format. * * iwsp(*) : (workspace) of declared length .ge. max(nrow,ncol). * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * Department of Mathematics, University of Queensland. * Brisbane QLD 4072, Australia. 1996. *----------------------------------------------------------------------| implicit none character from*3, to*3, diag*1 integer nrow, ncol, nz, ia(*), ja(*), iwsp(*) double precision a(*) integer i, j, k, nn, iflag character cpy_from*3, cpy_to*3, cpy_diag*1, c*1, S1*3, S2*3 logical LSAME intrinsic CHAR,ICHAR,LLE,LGE,MIN LSAME(S1,S2) = LLE(S1,S2).and.LGE(S1,S2) * *--- upper case strings ... * do k = 1,3 c = from(k:k) if ( c.ge.'a' .and. c.le.'z' ) c = CHAR( ICHAR(c)-32 ) cpy_from(k:k) = c c = to(k:k) if ( c.ge.'a' .and. c.le.'z' ) c = CHAR( ICHAR(c)-32 ) cpy_to(k:k) = c end do c = diag(1:1) if ( c.ge.'a' .and. c.le.'z' ) c = CHAR( ICHAR(c)-32 ) cpy_diag = c *--- quick return if possible ... if ( LSAME(cpy_from,cpy_to) .and. cpy_diag.eq.'N' ) return iflag = 1 if ( LSAME(cpy_from,'COO') ) iflag = 0 if ( LSAME(cpy_from,'CRS') ) iflag = 0 if ( LSAME(cpy_from,'CCS') ) iflag = 0 if ( iflag.ne.0 ) stop 'unexpected i/o formats (in DGCNVR)' iflag = 1 if ( LSAME(cpy_to,'COO') ) iflag = 0 if ( LSAME(cpy_to,'CRS') ) iflag = 0 if ( LSAME(cpy_to,'CCS') ) iflag = 0 if ( iflag.ne.0 ) stop 'unexpected i/o formats (in DGCNVR)' * *--- transit via COOrdinates format if input is not in COO ... * if ( LSAME(cpy_from,'CRS') ) then *--- expand ia indices ... nn = nz do i = nrow,1,-1 do k = 1,ia(i+1)-ia(i) ia(nn) = i nn = nn-1 end do end do endif if ( LSAME(cpy_from,'CCS') ) then *--- expand ja indices ... nn = nz do j = ncol,1,-1 do k = 1,ja(j+1)-ja(j) ja(nn) = j nn = nn-1 end do end do endif * *-- if requested, insert diagonal elements even if they are zero... * if ( cpy_diag.eq.'D' ) then nn = MIN( nrow,ncol ) do i = 1,nn iwsp(i) = 0 end do do k = 1,nz if ( ia(k).eq.ja(k) ) iwsp(ia(k)) = 1 end do do i = 1,nn if ( iwsp(i).eq.0 ) then nz = nz + 1 ia(nz) = i ja(nz) = i a(nz) = 0.0d0 endif end do endif *--- COO convertion ... if ( LSAME(cpy_to,'COO') ) return *--- CRS convertion ... if ( LSAME(cpy_to,'CRS') ) call dcmpac( nrow,nz,ia,ja,a, iwsp ) *--- CCS convertion ... if ( LSAME(cpy_to,'CCS') ) call dcmpac( ncol,nz,ja,ia,a, iwsp ) END subroutine dgcoov ( x, y ) c*********************************************************************72 c cc DGCOOV computes y = A * x, with A stored in "COO" format. c * *--- Computes y = A*x. A is passed via a fortran `common statement'. *--- A is assumed here to be under the COOrdinates storage format. * *----------------------------------------------------------------------| *-------------------------------NOTE-----------------------------------* * This is an accessory to Expokit and it is not intended to be * * complete. It is supplied primarily to ensure an unconstrained * * distribution and portability of the package. The matrix-vector * * multiplication routines supplied here fit the non symmetric * * storage and for a symmetric matrix, the entire (not half) matrix * * is required. If the sparsity pattern is known a priori, it is * * recommended to use the most advantageous format and to devise * * the most advantageous matrix-vector multiplication routine. * *----------------------------------------------------------------------* *----------------------------------------------------------------------* implicit none double precision x(*), y(*) integer n, nz, nzmax parameter( nzmax = 600000 ) integer ia(nzmax), ja(nzmax) double precision a(nzmax) common /RMAT/ a, ia, ja, nz, n integer i, j do j = 1,n y(j) = 0.0d0 end do do i = 1,nz y(ia(i)) = y(ia(i)) + a(i)*x(ja(i)) end do END subroutine dgcrsv ( x, y ) c*********************************************************************72 c cc DGCRSV computes y = A * x, for A stored in Compressed Row Storage (CRS) format. c implicit none double precision x(*), y(*) * *--- Computes y = A*x. A is passed via a fortran `common statement'. *--- A is assumed to be under the Compress Row Storage (CRS) format. * integer n, nz, nzmax parameter( nzmax = 600000 ) integer ia(nzmax), ja(nzmax) double precision a(nzmax) common /RMAT/ a, ia, ja, nz, n integer i, j do i = 1,n y(i) = 0.0d0 do j = ia(i),ia(i+1)-1 y(i) = y(i) + a(j)*x(ja(j)) end do end do END subroutine dgefa ( a, lda, n, ipvt, info ) c*********************************************************************72 c cc DGEFA factors a double precision matrix by gaussian elimination. c implicit none integer lda,n,ipvt(n),info double precision a(lda,n) c c dgefa is usually called by dgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for dgeco) = (1 + 9/n)*(time for dgefa) . c c on entry c c a double precision(lda, n) c the matrix to be factored. 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 and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this routine, but it does c indicate that dgesl or dgedi will divide by zero c if called. use rcond in dgeco for a reliable c indication of singularity. 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 daxpy,dscal,idamax c c internal variables c double precision t integer idamax,j,k,kp1,l,nm1 c c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = idamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0d0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0d0/a(k,k) call dscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0d0) info = n return end subroutine dgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, & beta, c, ldc ) c*********************************************************************72 c cc DGEMM performs matrix-matrix operations of the form C=alpha*A*B+beta*C. c implicit none CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * DGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN END subroutine dgemv ( trans, m, n, alpha, a, lda, x, incx, & beta, y, incy ) c*********************************************************************72 c cc DGEMV performs a matrix-vector operation of the form y=alpha*A*x+beta*y. c implicit none DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of DGEMV . * END subroutine dgesl ( a, lda, n, ipvt, b, job ) c*********************************************************************72 c cc DGESL solves a factored linear system of the form A*x=b. c implicit none integer lda,n,ipvt(n),job double precision a(lda,n),b(n) c c dgesl solves the double precision system c a * x = b or trans(a) * x = b c using the factors computed by dgeco or dgefa. c c on entry c c a double precision(lda, n) c the output from dgeco or dgefa. 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 ipvt integer(n) c the pivot vector from dgeco or dgefa. c c b double precision(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if dgeco has set rcond .gt. 0.0 c or dgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call dgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call dgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue 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 daxpy,ddot c c internal variables c double precision ddot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call daxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call daxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = ddot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine dgesv ( n, m, a,lda, ipiv, b,ldb, iflag ) c*********************************************************************72 c cc DGESV solves multiple general linear systems of the form A*X=B. c * This is a lightweight substitute to the external LAPACK routines * used by EXPOKIT. It is supplied to ensure that EXPOKIT is * self-contained and can still run if LAPACK is not yet installed * in your environement. c implicit none integer N, M, LDA, LDB, IPIV(N), IFLAG double precision A(LDA,N), B(LDB,M) integer j call DGEFA( A,LDA, N, IPIV, IFLAG ) if ( IFLAG.ne.0 ) stop "Error in DGESV (LU factorisation)" do j = 1,M call DGESL( A,LDA, N, IPIV,B(1,j), 0 ) end do end subroutine dgexpv ( n, m, t, v, w, tol, anorm, . wsp, lwsp, iwsp, liwsp, matvec, itrace, iflag ) c*********************************************************************72 c cc DGEXPV computes w=exp(t*A)*v for a general matrix A. c implicit none integer n, m, lwsp, liwsp, itrace, iflag, iwsp(liwsp) double precision t, tol, anorm, v(n), w(n), wsp(lwsp) external matvec *-----Purpose----------------------------------------------------------| * *--- DGEXPV computes w = exp(t*A)*v - for a General matrix A. * * It does not compute the matrix exponential in isolation but * instead, it computes directly the action of the exponential * operator on the operand vector. This way of doing so allows * for addressing large sparse problems. * * The method used is based on Krylov subspace projection * techniques and the matrix under consideration interacts only * via the external routine `matvec' performing the matrix-vector * product (matrix-free method). * *-----Arguments--------------------------------------------------------| * * n : (input) order of the principal matrix A. * * m : (input) maximum size for the Krylov basis. * * t : (input) time at wich the solution is needed (can be < 0). * * v(n) : (input) given operand vector. * * w(n) : (output) computed approximation of exp(t*A)*v. * * tol : (input/output) the requested accuracy tolerance on w. * If on input tol=0.0d0 or tol is too small (tol.le.eps) * the internal value sqrt(eps) is used, and tol is set to * sqrt(eps) on output (`eps' denotes the machine epsilon). * (`Happy breakdown' is assumed if h(j+1,j) .le. anorm*tol) * * anorm : (input) an approximation of some norm of A. * * wsp(lwsp): (workspace) lwsp .ge. n*(m+1)+n+(m+2)^2+4*(m+2)^2+ideg+1 * +---------+-------+---------------+ * (actually, ideg=6) V H wsp for PADE * * iwsp(liwsp): (workspace) liwsp .ge. m+2 * * matvec : external routine for matrix-vector multiplication. * synopsis: matvec( x, y ) * double precision x(*), y(*) * computes: y(1:n) <- A*x(1:n) * where A is the principal matrix. * * itrace : (input) running mode. 0=silent, 1=print step-by-step info * * iflag : (output) exit flag. * <0 - bad input arguments * 0 - no problem * 1 - maximum number of steps reached without convergence * 2 - requested tolerance was too high * *-----Accounts on the computation--------------------------------------| * Upon exit, an interested user may retrieve accounts on the * computations. They are located in wsp and iwsp as indicated below: * * location mnemonic description * -----------------------------------------------------------------| * iwsp(1) = nmult, number of matrix-vector multiplications used * iwsp(2) = nexph, number of Hessenberg matrix exponential evaluated * iwsp(3) = nscale, number of repeated squaring involved in Pade * iwsp(4) = nstep, number of integration steps used up to completion * iwsp(5) = nreject, number of rejected step-sizes * iwsp(6) = ibrkflag, set to 1 if `happy breakdown' and 0 otherwise * iwsp(7) = mbrkdwn, if `happy brkdown', basis-size when it occured * -----------------------------------------------------------------| * wsp(1) = step_min, minimum step-size used during integration * wsp(2) = step_max, maximum step-size used during integration * wsp(3) = dummy * wsp(4) = dummy * wsp(5) = x_error, maximum among all local truncation errors * wsp(6) = s_error, global sum of local truncation errors * wsp(7) = tbrkdwn, if `happy breakdown', time when it occured * wsp(8) = t_now, integration domain successfully covered * wsp(9) = hump, i.e., max||exp(sA)||, s in [0,t] (or [t,0] if t<0) * wsp(10) = ||w||/||v||, scaled norm of the solution w. * -----------------------------------------------------------------| * The `hump' is a measure of the conditioning of the problem. The * matrix exponential is well-conditioned if hump = 1, whereas it is * poorly-conditioned if hump >> 1. However the solution can still be * relatively fairly accurate even when the hump is large (the hump * is an upper bound), especially when the hump and the scaled norm * of w [this is also computed and returned in wsp(10)] are of the * same order of magnitude (further details in reference below). * *----------------------------------------------------------------------| *-----The following parameters may also be adjusted herein-------------| * integer mxstep, mxreject, ideg double precision delta, gamma parameter( mxstep = 1000, . mxreject = 0, . ideg = 6, . delta = 1.2d0, . gamma = 0.9d0 ) * mxstep : maximum allowable number of integration steps. * The value 0 means an infinite number of steps. * * mxreject: maximum allowable number of rejections at each step. * The value 0 means an infinite number of rejections. * * ideg : the Pade approximation of type (ideg,ideg) is used as * an approximation to exp(H). The value 0 switches to the * uniform rational Chebyshev approximation of type (14,14) * * delta : local truncation error `safety factor' * * gamma : stepsize `shrinking factor' * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * EXPOKIT: Software Package for Computing Matrix Exponentials. * ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 *----------------------------------------------------------------------| * integer i, j, k1, mh, mx, iv, ih, j1v, ns, ifree, lfree, iexph, . ireject,ibrkflag,mbrkdwn, nmult, nreject, nexph, nscale, . nstep double precision sgn, t_out, tbrkdwn, step_min,step_max, err_loc, . s_error, x_error, t_now, t_new, t_step, t_old, . xm, beta, break_tol, p1, p2, p3, eps, rndoff, . vnorm, avnorm, hj1j, hij, hump, SQR1 intrinsic AINT,ABS,DBLE,LOG10,MAX,MIN,NINT,SIGN,SQRT double precision DDOT, DNRM2 *--- check restrictions on input parameters ... iflag = 0 if ( lwsp.lt.n*(m+2)+5*(m+2)**2+ideg+1 ) iflag = -1 if ( liwsp.lt.m+2 ) iflag = -2 if ( m.ge.n .or. m.le.0 ) iflag = -3 if ( iflag.ne.0 ) stop 'bad sizes (in input of DGEXPV)' * *--- initialisations ... * k1 = 2 mh = m + 2 iv = 1 ih = iv + n*(m+1) + n ifree = ih + mh*mh lfree = lwsp - ifree + 1 ibrkflag = 0 mbrkdwn = m nmult = 0 nreject = 0 nexph = 0 nscale = 0 t_out = ABS( t ) tbrkdwn = 0.0d0 step_min = t_out step_max = 0.0d0 nstep = 0 s_error = 0.0d0 x_error = 0.0d0 t_now = 0.0d0 t_new = 0.0d0 p1 = 4.0d0/3.0d0 1 p2 = p1 - 1.0d0 p3 = p2 + p2 + p2 eps = ABS( p3-1.0d0 ) if ( eps.eq.0.0d0 ) go to 1 if ( tol.le.eps ) tol = SQRT( eps ) rndoff = eps*anorm break_tol = 1.0d-7 *>>> break_tol = tol *>>> break_tol = anorm*tol sgn = SIGN( 1.0d0,t ) call DCOPY( n, v,1, w,1 ) beta = DNRM2( n, w,1 ) vnorm = beta hump = beta * *--- obtain the very first stepsize ... * SQR1 = SQRT( 0.1d0 ) xm = 1.0d0/DBLE( m ) p1 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p1/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 * *--- step-by-step integration ... * 100 if ( t_now.ge.t_out ) goto 500 nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) p1 = 1.0d0/beta do i = 1,n wsp(iv + i-1) = p1*w(i) end do do i = 1,mh*mh wsp(ih+i-1) = 0.0d0 end do * *--- Arnoldi loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) do i = 1,j hij = DDOT( n, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) call DAXPY( n, -hij, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) wsp(ih+(j-1)*mh+i-1) = hij end do hj1j = DNRM2( n, wsp(j1v),1 ) *--- if `happy breakdown' go straightforward at the end ... if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = hj1j call DSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DNRM2( n, wsp(j1v),1 ) * *--- set 1 for the 2-corrected scheme ... * 300 continue wsp(ih+m*mh+m+1) = 1.0d0 * *--- loop while ireject>> break_tol = tol *>>> break_tol = anorm*tol * *--- step-by-step integration ... * sgn = SIGN( 1.0d0,t ) SQR1 = SQRT( 0.1d0 ) call DCOPY( n, v,1, w,1 ) 100 if ( t_now.ge.t_out ) goto 500 nmult = nmult + 1 call matvec( w, wsp(iv) ) call DAXPY( n, 1.0d0, u,1, wsp(iv),1 ) beta = DNRM2( n, wsp(iv),1 ) if ( beta.eq.0.0d0 ) goto 500 call DSCAL( n, 1.0d0/beta, wsp(iv),1 ) do i = 1,mh*mh wsp(ih+i-1) = 0.0d0 end do if ( nstep.eq.0 ) then *--- obtain the very first stepsize ... xm = 1.0d0/DBLE( m ) p1 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p1/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 endif nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) * *--- Arnoldi loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) do i = 1,j hij = DDOT( n, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) call DAXPY( n, -hij, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) wsp(ih+(j-1)*mh+i-1) = hij end do hj1j = DNRM2( n, wsp(j1v),1 ) *--- if `happy breakdown' go straightforward at the end ... if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = hj1j call DSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DNRM2( n, wsp(j1v),1 ) * *--- set 1's for the 3-extended scheme ... * 300 continue wsp(ih+mh*mbrkdwn) = 1.0d0 wsp(ih+(m-1)*mh+m) = 0.0d0 do i = 1,k1-1 wsp(ih+(m+i)*mh+m+i-1) = 1.0d0 end do * *--- loop while ireject> 1. However the solution can still be * relatively fairly accurate even when the hump is large (the hump * is an upper bound), especially when the hump and the scaled norm * of w [this is also computed and returned in wsp(10)] are of the * same order of magnitude (further details in reference below). * Markov chains are usually well-conditioned problems. * *----------------------------------------------------------------------| *-----The following parameters may also be adjusted herein-------------| * integer mxstep, mxreject, ideg double precision delta, gamma parameter( mxstep = 500, . mxreject = 0, . ideg = 6, . delta = 1.2d0, . gamma = 0.9d0 ) * mxstep : maximum allowable number of integration steps. * The value 0 means an infinite number of steps. * * mxreject: maximum allowable number of rejections at each step. * The value 0 means an infinite number of rejections. * * ideg : the Pade approximation of type (ideg,ideg) is used as * an approximation to exp(H). The value 0 switches to the * uniform rational Chebyshev approximation of type (14,14) * * delta : local truncation error `safety factor' * * gamma : stepsize `shrinking factor' * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * EXPOKIT: Software Package for Computing Matrix Exponentials. * ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 *----------------------------------------------------------------------| * integer i, j, k1, mh, mx, iv, ih, j1v, ns, ifree, lfree, iexph, . ireject,ibrkflag,mbrkdwn, nmult, nreject, nexph, nscale, . nstep double precision sgn, t_out, tbrkdwn, step_min,step_max, err_loc, . s_error, x_error, t_now, t_new, t_step, t_old, . xm, beta, break_tol, p1, p2, p3, eps, rndoff, . vnorm, avnorm, hj1j, hij, hump, SQR1, . roundoff, s_round, x_round intrinsic AINT,ABS,DBLE,LOG10,MAX,MIN,NINT,SIGN,SQRT double precision DDOT, DNRM2, DASUM *--- check restrictions on input parameters ... iflag = 0 if ( lwsp.lt.n*(m+2)+5*(m+2)**2+ideg+1 ) iflag = -1 if ( liwsp.lt.m+2 ) iflag = -2 if ( m.ge.n .or. m.le.0 ) iflag = -3 if ( iflag.ne.0 ) stop 'bad sizes (in input of DMEXPV)' * *--- initialisations ... * k1 = 2 mh = m + 2 iv = 1 ih = iv + n*(m+1) + n ifree = ih + mh*mh lfree = lwsp - ifree + 1 ibrkflag = 0 mbrkdwn = m nmult = 0 nreject = 0 nexph = 0 nscale = 0 sgn = SIGN( 1.0d0,t ) t_out = ABS( t ) tbrkdwn = 0.0d0 step_min = t_out step_max = 0.0d0 nstep = 0 s_error = 0.0d0 s_round = 0.0d0 x_error = 0.0d0 x_round = 0.0d0 t_now = 0.0d0 t_new = 0.0d0 p1 = 4.0d0/3.0d0 1 p2 = p1 - 1.0d0 p3 = p2 + p2 + p2 eps = ABS( p3-1.0d0 ) if ( eps.eq.0.0d0 ) go to 1 if ( tol.le.eps ) tol = SQRT( eps ) rndoff = eps*anorm break_tol = 1.0d-7 *>>> break_tol = tol *>>> break_tol = anorm*tol call DCOPY( n, v,1, w,1 ) beta = DNRM2( n, w,1 ) vnorm = beta hump = beta * *--- obtain the very first stepsize ... * SQR1 = SQRT( 0.1d0 ) xm = 1.0d0/DBLE( m ) p1 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p1/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 * *--- step-by-step integration ... * 100 if ( t_now.ge.t_out ) goto 500 nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) p1 = 1.0d0/beta do i = 1,n wsp(iv + i-1) = p1*w(i) end do do i = 1,mh*mh wsp(ih+i-1) = 0.0d0 end do * *--- Arnoldi loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) do i = 1,j hij = DDOT( n, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) call DAXPY( n, -hij, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) wsp(ih+(j-1)*mh+i-1) = hij end do hj1j = DNRM2( n, wsp(j1v),1 ) *--- if `happy breakdown' go straightforward at the end ... if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = hj1j call DSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DNRM2( n, wsp(j1v),1 ) * *--- set 1 for the 2-corrected scheme ... * 300 continue wsp(ih+m*mh+m+1) = 1.0d0 * *--- loop while ireject> 1. However the solution can still be * relatively fairly accurate even when the hump is large (the hump * is an upper bound), especially when the hump and the scaled norm * of w [this is also computed and returned in wsp(10)] are of the * same order of magnitude (further details in reference below). * *----------------------------------------------------------------------| *-----The following parameters may also be adjusted herein-------------| * integer mxstep, mxreject, ideg double precision delta, gamma parameter( mxstep = 500, . mxreject = 0, . ideg = 6, . delta = 1.2d0, . gamma = 0.9d0 ) * mxstep : maximum allowable number of integration steps. * The value 0 means an infinite number of steps. * * mxreject: maximum allowable number of rejections at each step. * The value 0 means an infinite number of rejections. * * ideg : the Pade approximation of type (ideg,ideg) is used as * an approximation to exp(H). The value 0 switches to the * uniform rational Chebyshev approximation of type (14,14) * * delta : local truncation error `safety factor' * * gamma : stepsize `shrinking factor' * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * EXPOKIT: Software Package for Computing Matrix Exponentials. * ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 *----------------------------------------------------------------------| * integer i, j, k1, mh, mx, iv, ih, j1v, ns, ifree, lfree, iexph, . ireject,ibrkflag,mbrkdwn, nmult, nreject, nexph, nscale, . nstep double precision sgn, t_out, tbrkdwn, step_min,step_max, err_loc, . s_error, x_error, t_now, t_new, t_step, t_old, . xm, beta, break_tol, p1, p2, p3, eps, rndoff, . vnorm, avnorm, hj1j, hjj, hump, SQR1 intrinsic AINT,ABS,DBLE,LOG10,MAX,MIN,NINT,SIGN,SQRT double precision DDOT, DNRM2 *--- check restrictions on input parameters ... iflag = 0 if ( lwsp.lt.n*(m+2)+5*(m+2)**2+ideg+1 ) iflag = -1 if ( liwsp.lt.m+2 ) iflag = -2 if ( m.ge.n .or. m.le.0 ) iflag = -3 if ( iflag.ne.0 ) stop 'bad sizes (in input of DSEXPV)' * *--- initialisations ... * k1 = 2 mh = m + 2 iv = 1 ih = iv + n*(m+1) + n ifree = ih + mh*mh lfree = lwsp - ifree + 1 ibrkflag = 0 mbrkdwn = m nmult = 0 nreject = 0 nexph = 0 nscale = 0 t_out = ABS( t ) tbrkdwn = 0.0d0 step_min = t_out step_max = 0.0d0 nstep = 0 s_error = 0.0d0 x_error = 0.0d0 t_now = 0.0d0 t_new = 0.0d0 p1 = 4.0d0/3.0d0 1 p2 = p1 - 1.0d0 p3 = p2 + p2 + p2 eps = ABS( p3-1.0d0 ) if ( eps.eq.0.0d0 ) go to 1 if ( tol.le.eps ) tol = SQRT( eps ) rndoff = eps*anorm break_tol = 1.0d-7 *>>> break_tol = tol *>>> break_tol = anorm*tol sgn = SIGN( 1.0d0,t ) call DCOPY( n, v,1, w,1 ) beta = DNRM2( n, w,1 ) vnorm = beta hump = beta * *--- obtain the very first stepsize ... * SQR1 = SQRT( 0.1d0 ) xm = 1.0d0/DBLE( m ) p1 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p1/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 * *--- step-by-step integration ... * 100 if ( t_now.ge.t_out ) goto 500 nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) p1 = 1.0d0/beta do i = 1,n wsp(iv + i-1) = p1*w(i) end do do i = 1,mh*mh wsp(ih+i-1) = 0.0d0 end do * *--- Lanczos loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) if ( j.gt.1 ) . call DAXPY(n,-wsp(ih+(j-1)*mh+j-2),wsp(j1v-2*n),1,wsp(j1v),1) hjj = DDOT( n, wsp(j1v-n),1, wsp(j1v),1 ) call DAXPY( n, -hjj, wsp(j1v-n),1, wsp(j1v),1 ) hj1j = DNRM2( n, wsp(j1v),1 ) wsp(ih+(j-1)*(mh+1)) = hjj *--- if `happy breakdown' go straightforward at the end ... if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = hj1j wsp(ih+j*mh+j-1) = hj1j call DSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DNRM2( n, wsp(j1v),1 ) * *--- set 1 for the 2-corrected scheme ... * 300 continue wsp(ih+m*mh+m-1) = 0.0d0 wsp(ih+m*mh+m+1) = 1.0d0 * *--- loop while ireject>> break_tol = tol *>>> break_tol = anorm*tol * *--- step-by-step integration ... * sgn = SIGN( 1.0d0,t ) SQR1 = SQRT( 0.1d0 ) call DCOPY( n, v,1, w,1 ) 100 if ( t_now.ge.t_out ) goto 500 nmult = nmult + 1 call matvec( w, wsp(iv) ) call DAXPY( n, 1.0d0, u,1, wsp(iv),1 ) beta = DNRM2( n, wsp(iv),1 ) if ( beta.eq.0.0d0 ) goto 500 call DSCAL( n, 1.0d0/beta, wsp(iv),1 ) do i = 1,mh*mh wsp(ih+i-1) = 0.0d0 end do if ( nstep.eq.0 ) then *--- obtain the very first stepsize ... xm = 1.0d0/DBLE( m ) p1 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p1/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 endif nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) * *--- Lanczos loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) if ( j.gt.1 ) . call DAXPY(n,-wsp(ih+(j-1)*mh+j-2),wsp(j1v-2*n),1,wsp(j1v),1) hjj = DDOT( n, wsp(j1v-n),1, wsp(j1v),1 ) call DAXPY( n, -hjj, wsp(j1v-n),1, wsp(j1v),1 ) hj1j = DNRM2( n, wsp(j1v),1 ) wsp(ih+(j-1)*(mh+1)) = hjj *--- if `happy breakdown' go straightforward at the end ... if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = hj1j wsp(ih+j*mh+j-1) = hj1j call DSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DNRM2( n, wsp(j1v),1 ) * *--- set 1's for the 3-extended scheme ... * 300 continue wsp(ih+mh*mbrkdwn) = 1.0d0 wsp(ih+m*mh+m-1) = 0.0d0 wsp(ih+(m-1)*mh+m) = 0.0d0 do i = 1,k1-1 wsp(ih+(m+i)*mh+m+i-1) = 1.0d0 end do * *--- loop while ireject= nz. If the matrix is symmetric, both * lower and upper parts are included explicitly * * iwsp (workspace) of length >= n * *----------------------------------------------------------------------| * character title*72, key*8, type*3, ptrfmt*16, . indfmt*16, valfmt*20, rhsfmt*20, rhstyp*1 integer totcrd, ptrcrd, indcrd, valcrd, rhscrd, nrow, . nrhs, nrhsix, i, j, k, io, nn, nnz intrinsic INDEX open( UNIT=7, STATUS='old', IOstat=io, FILE=filename ) if ( io.ne.0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LOADHB - Fatal error!' write ( *, '(a)' ) ' Could not access Harwell-Boeing matrix' write ( *, '(a)' ) ' "' // trim ( filename ) // '".' stop end if read( UNIT=7,FMT=10 ) title, key, . totcrd, ptrcrd, indcrd, valcrd, rhscrd, . type, nrow, nn, nnz, nrhs, . ptrfmt, indfmt, valfmt, rhsfmt print*, title,'type :',type,' size :',nrow,nn print*,'order :',nn,' number of nonzero :',nnz if ( nn.gt.n ) stop 'in LOADHB. Please increase n' if ( nnz.gt.nz ) stop 'in LOADHB. Please increase nz' *--- leave if there is no values ... if ( valcrd.le.0 ) stop 'Empty Harwell-Boeing matrix' if ( rhscrd.gt.0 ) then read( UNIT=7,FMT=11 ) rhstyp, nrhs, nrhsix print*,'There is a second hand' endif *--- read data... read( UNIT=7,FMT=ptrfmt ) (ja(i), i = 1,nn+1) read( UNIT=7,FMT=indfmt ) (ia(i), i = 1,nnz) read( UNIT=7,FMT=valfmt ) (a(i), i = 1,nnz) close( UNIT=7 ) *--- for the sake of experiments, store both parts if symmetric matrix if ( type.eq.'RSA' ) then *--- expand ja indices ... k = nnz do j = nn,1,-1 do i = 1,ja(j+1)-ja(j) ja(k) = j k = k-1 end do end do *--- insert the other half ... k = nnz do i = 1,k if ( ia(i).ne.ja(i) ) then nnz = nnz + 1 if ( nnz.gt.nz ) stop 'in LOADHB. Please increase nz' ia(nnz) = ja(i) ja(nnz) = ia(i) a(nnz) = a(i) endif end do type = 'COO' else type = 'CCS' endif call dgcnvr( type,spformat,'n', nn,nn, nnz, ia, ja, a, iwsp ) n = nn nz = nnz print*,'Harwell-Boeing matrix loaded' 10 format(A72, A8/ 5i14 / A3, 11x, 4i14 / 2a16, 2a20) 11 format(A1, 13x, 2i14) end function lsame( ca, cb ) c*********************************************************************72 c cc LSAME compares two characters for identity. * * -- LAPACK auxiliary routine (version 1.1) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * implicit none logical lsame CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Discussion: c c This FORTRAN77 version is made available for cases where the c FORTRAN90 version cannot be used. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, month(m), y, h, ':', n, ':', s, '.', mm, ampm return end subroutine tnspos ( n, nz, ia, ja, a, iwsp ) c*********************************************************************72 c cc TNSPOS transposes a CRS matrix. c implicit none integer n, nz, ia(nz), ja(nz), iwsp(n) double precision a(nz) integer i, j call dgcnvr( 'crs','ccs','n', n,n, nz, ia, ja, a, iwsp ) do i = 1,nz j = ja(i) ja(i) = ia(i) ia(i) = j end do END subroutine xerbla ( srname, info ) c*********************************************************************72 c cc XERBLA is an error handler for the level 2 BLAS routines. c implicit none INTEGER INFO CHARACTER*6 SRNAME * .. * * Purpose * ======= * * XERBLA is an error handler for the Level 2 BLAS routines. * * It is called by the Level 2 BLAS routines if an input parameter is * invalid. * * Installers should consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Parameters * ========== * * SRNAME - CHARACTER*6. * On entry, SRNAME specifies the name of the routine which * called XERBLA. * * INFO - INTEGER. * On entry, INFO specifies the position of the invalid * parameter in the parameter-list of the calling routine. * * * Auxiliary routine for Level 2 Blas. * * Written on 20-July-1986. * * .. Executable Statements .. * WRITE (*,99999) SRNAME, INFO * STOP * 99999 FORMAT ( ' ** On entry to ', A6, ' parameter number ', I2, $ ' had an illegal value' ) * * End of XERBLA. * END subroutine zaxpy ( n, za, zx, incx, zy, incy ) c*********************************************************************72 c cc ZAXPY adds a multiple of one vector to another. c c jack dongarra, 3/11/78. c implicit none double precision dcabs1 integer i integer incx integer incy integer ix integer iy integer n complex*16 za complex*16 zx(1) complex*16 zy(1) if(n.le.0)return if (dcabs1(za) .eq. 0.0d0) return if (incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n zy(iy) = zy(iy) + za*zx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n zy(i) = zy(i) + za*zx(i) 30 continue return end subroutine zcmpac ( n, nx, ix, ixx, xx, iwsp ) c*********************************************************************72 c cc ZCMPAC compacts the array IX and sorts IXX and XX. c *-- (This is a gateway routine for ZGCNVR) ... c implicit none integer n, nx, ix(nx), ixx(nx), iwsp(n) complex*16 xx(nx) integer k * *--- sort ix and carry ixx and xx along ... * call izsrt2( nx, ix, ixx, xx ) * *--- adjust pointers ... * do k = 1,n iwsp(k) = 0 end do do k = 1,nx iwsp(ix(k)) = iwsp(ix(k)) + 1 end do ix(n+1) = nx + 1 do k = n,1,-1 ix(k) = ix(k+1)-iwsp(k) end do * *--- sort ixx in increasing order and carry xx along ... * do k = 1,n call izsrt1( iwsp(k), ixx(ix(k)), xx(ix(k)) ) end do END subroutine zcopy ( n, zx, incx, zy, incy ) c*********************************************************************72 c cc ZCOPY copies a vector. c c jack dongarra, linpack, 4/11/78. c implicit none double complex zx(1),zy(1) integer i,incx,incy,ix,iy,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n zy(iy) = zx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n zy(i) = zx(i) 30 continue return end function zdotc(n,zx,incx,zy,incy) c*********************************************************************72 c cc ZDOTC computes the conjugated dot product of complex vectors. c c jack dongarra, 3/11/78. c implicit none integer i integer incx integer incy integer ix integer iy integer n double complex zdotc double complex zx(1),zy(1),ztemp ztemp = (0.0d0,0.0d0) zdotc = (0.0d0,0.0d0) if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n ztemp = ztemp + dconjg(zx(ix))*zy(iy) ix = ix + incx iy = iy + incy 10 continue zdotc = ztemp return c c code for both increments equal to 1 c 20 do 30 i = 1,n ztemp = ztemp + dconjg(zx(i))*zy(i) 30 continue zdotc = ztemp return end function zdotu(n,zx,incx,zy,incy) c*********************************************************************72 c cc ZDOTU computes the unconjugated dot product of complex vectors. c c jack dongarra, 3/11/78. c implicit none integer i integer incx integer incy integer ix integer iy integer n double complex zdotu double complex zx(1),zy(1),ztemp ztemp = (0.0d0,0.0d0) zdotu = (0.0d0,0.0d0) if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n ztemp = ztemp + zx(ix)*zy(iy) ix = ix + incx iy = iy + incy 10 continue zdotu = ztemp return c c code for both increments equal to 1 c 20 do 30 i = 1,n ztemp = ztemp + zx(i)*zy(i) 30 continue zdotu = ztemp return end subroutine zdscal ( n, da, zx, incx ) c*********************************************************************72 c cc ZDSCAL scales a complex vector by a real constant. c c jack dongarra, 3/11/78. c modified 3/93 to return if incx .le. 0. c implicit none double complex zx(1) double precision da integer i,incx,ix,n c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 do 10 i = 1,n zx(ix) = dcmplx(da,0.0d0)*zx(ix) ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 do 30 i = 1,n zx(i) = dcmplx(da,0.0d0)*zx(i) 30 continue return end subroutine zgccsv ( x, y ) c*********************************************************************72 c cc ZGCCSV computes y = A * x, for A in compressed column storage (CCS). c implicit none complex*16 x(*), y(*) * *--- Computes y = A*x. A is passed via a fortran `common statement'. *--- A is assumed to be under the Compress Column Storage (CCS) format. * integer n, nz, nzmax parameter( nzmax = 50000 ) integer ia(nzmax), ja(nzmax) complex*16 a(nzmax) common /CMAT/ a, ia, ja, nz, n integer i, j complex*16 ZERO parameter( ZERO=(0.0d0,0.0d0) ) do i = 1,n y(i) = ZERO end do do j = 1,n do i = ja(j),ja(j+1)-1 y(ia(i)) = y(ia(i)) + a(i)*x(j) end do end do end subroutine zgchbv ( m, t, h,ldh, y, wsp, iwsp, iflag ) c*********************************************************************72 c cc ZGCHBV computes y = exp(t*H)*y using a partial fraction expansion. c implicit none integer m, ldh, iflag, iwsp(m) double precision t complex*16 H(ldh,m), y(m), wsp(m*(m+2)) *-----Purpose----------------------------------------------------------| * *--- ZGCHBV computes y = exp(t*H)*y using the partial fraction * expansion of the uniform rational Chebyshev approximation * to exp(-x) of type (14,14). H is a General matrix. * About 14-digit accuracy is expected if the matrix H is negative * definite. The algorithm may behave poorly otherwise. * *-----Arguments--------------------------------------------------------| * * m : (input) order of the matrix H. * * t : (input) time-scaling factor (can be < 0). * * H(ldh,m): (input) argument matrix. * * y(m) : (input/output) on input the operand vector, * on output the resulting vector exp(t*H)*y. * * iwsp(m) : (workspace) * * wsp : (workspace). Observe that a double precision vector of * length 2*m*(m+2) can be used as well when calling this * routine. * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * EXPOKIT: Software Package for Computing Matrix Exponentials. * ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 *----------------------------------------------------------------------| * integer ndeg, i, j, ip, ih, iy, iz parameter ( ndeg=7 ) double precision alpha0 complex*16 alpha(2*ndeg), theta(2*ndeg) *--- Pointers ... ih = 1 iy = ih + m*m iz = iy + m *--- Coefficients and poles of the partial fraction expansion ... alpha0 = 0.183216998528140087D-11 alpha(1)=( 0.557503973136501826D+02,-0.204295038779771857D+03) alpha(2)=(-0.938666838877006739D+02, 0.912874896775456363D+02) alpha(3)=( 0.469965415550370835D+02,-0.116167609985818103D+02) alpha(4)=(-0.961424200626061065D+01,-0.264195613880262669D+01) alpha(5)=( 0.752722063978321642D+00, 0.670367365566377770D+00) alpha(6)=(-0.188781253158648576D-01,-0.343696176445802414D-01) alpha(7)=( 0.143086431411801849D-03, 0.287221133228814096D-03) theta(1)=(-0.562314417475317895D+01, 0.119406921611247440D+01) theta(2)=(-0.508934679728216110D+01, 0.358882439228376881D+01) theta(3)=(-0.399337136365302569D+01, 0.600483209099604664D+01) theta(4)=(-0.226978543095856366D+01, 0.846173881758693369D+01) theta(5)=( 0.208756929753827868D+00, 0.109912615662209418D+02) theta(6)=( 0.370327340957595652D+01, 0.136563731924991884D+02) theta(7)=( 0.889777151877331107D+01, 0.166309842834712071D+02) * do ip = 1,ndeg theta(ndeg+ip) = CONJG( theta(ip) ) alpha(ndeg+ip) = CONJG( alpha(ip) ) end do * *--- Accumulation of the contribution of each pole ... * do j = 1,m wsp(iz+j-1) = y(j) y(j) = y(j)*alpha0 end do do ip = 1,2*ndeg alpha(ip) = 0.5d0*alpha(ip) *--- Solve each fraction using Gaussian elimination with pivoting... do j = 1,m do i = 1,m wsp(ih+(j-1)*m+i-1) = -t*H(i,j) end do wsp(ih+(j-1)*m+j-1) = wsp(ih+(j-1)*m+j-1)-theta(ip) wsp(iy+j-1) = wsp(iz+j-1) end do call ZGESV( M, 1, WSP(iH),M, IWSP, WSP(iY),M, IFLAG ) if ( IFLAG.ne.0 ) stop 'Error in ZGCHBV' *--- Accumulate the partial result in y ... do i = 1,m y(i) = y(i) + alpha(ip)*wsp(iy+i-1) end do end do END subroutine zgcnvr ( from,to,diag, nrow,ncol, nz, ia, ja, a, iwsp ) c*********************************************************************72 c cc ZGCNVR transforms a sparse storage format into another sparse storage format. c implicit none character from*3, to*3, diag*1 integer nrow, ncol, nz, ia(*), ja(*), iwsp(*) complex*16 a(*) *-----Purpose----------------------------------------------------------| * *--- ZGCNVR transforms a sparse storage format into another sparse * storage format. The matrix can be rectangular. * This is the complex counterpart of DGCNVR. * *-----Arguments--------------------------------------------------------| * * from : (input, character*3) the storage format holding the * matrix on input. Accepted values are: * 'COO' : COOrdinates * 'CRS' : Compressed Row Storage * 'CCS' : Compressed Column Storage. * * to : (input, character*3) the storage format holding the * matrix on output. Same accepted values as above. * * diag : (input, character*1) specifies whether or not the * entire diagonal should be lodged, i.e., including * null diagonal entries. (This may be required by * some applications that make use of the locations * of diagonal elements.) * diag = 'N', no specific attention is given to diagonal * elements. * diag = 'D', the entire diagonal is lodged, including * null elements on the diagonal. * if from=to & diag='D', null diagonal entries are * explicitly inserted in the actual matrix. * * nrow : (input) number of rows in the matrix. * * ncol : (input) number of columns in the matrix. * * nz : (input/output) number of non-zeros elements. * If diag='D' and null diagonal entries are inserted, then * nz is updated on exit and contains the effective number * of entries stored. In what follows, nz' (read nz prime) * denotes the updated value of nz, nz <= nz' <= nz+n. * If diag='N' then nz'=nz. * * ia(*) : (input/output) of declared length .ge. nz'. * On input, * if from='CRS', ia(1:nrow+1) contains pointers for the * beginning of each row. * if from='COO', or 'CCS', ia(1:nz) contains row indices. * On output, * if to='CRS', ia(1:nrow+1) contains pointers for the * beginning of each row. * if to='COO' or 'CCS', ia(1:nz') contains row indices * in increasing order in each column. * * ja(*) : (input/output) of declared length .ge. nz'. * On input, * if from='CRS', ja(1:ncol+1) contains pointers for the * beginning of each column. * if from='COO' or 'CCS', ja(1:nz) contains col. indices. * On output, * if to='CRS', ja(1:ncol+1) contains pointers for the * beginning of each column. * if to='COO' or 'CCS', ja(1:nz') contains column indices * in increasing order in each row. * * a(*) : (input/output) On input, a(1:nz) fits the input format * and on output, a(1:nz') fits the output format. * * iwsp(*) : (workspace) of declared length .ge. max(nrow,ncol). * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * Department of Mathematics, University of Queensland. * Brisbane QLD 4072, Australia. 1996. *----------------------------------------------------------------------| integer i, j, k, nn, iflag character cpy_from*3, cpy_to*3, cpy_diag*1, c*1, S1*3, S2*3 logical LSAME intrinsic CHAR,ICHAR,LLE,LGE,MIN LSAME(S1,S2) = LLE(S1,S2).and.LGE(S1,S2) * *--- upper case strings ... * do k = 1,3 c = from(k:k) if ( c.ge.'a' .and. c.le.'z' ) c = CHAR( ICHAR(c)-32 ) cpy_from(k:k) = c c = to(k:k) if ( c.ge.'a' .and. c.le.'z' ) c = CHAR( ICHAR(c)-32 ) cpy_to(k:k) = c end do c = diag(1:1) if ( c.ge.'a' .and. c.le.'z' ) c = CHAR( ICHAR(c)-32 ) cpy_diag = c *--- quick return if possible ... if ( LSAME(cpy_from,cpy_to) .and. cpy_diag.eq.'N' ) return iflag = 1 if ( LSAME(cpy_from,'COO') ) iflag = 0 if ( LSAME(cpy_from,'CRS') ) iflag = 0 if ( LSAME(cpy_from,'CCS') ) iflag = 0 if ( iflag.ne.0 ) stop 'unexpected i/o formats (in ZGCNVR)' iflag = 1 if ( LSAME(cpy_to,'COO') ) iflag = 0 if ( LSAME(cpy_to,'CRS') ) iflag = 0 if ( LSAME(cpy_to,'CCS') ) iflag = 0 if ( iflag.ne.0 ) stop 'unexpected i/o formats (in ZGCNVR)' * *--- transit via COOrdinates format if input is not in COO ... * if ( LSAME(cpy_from,'CRS') ) then *--- expand ia indices ... nn = nz do i = nrow,1,-1 do k = 1,ia(i+1)-ia(i) ia(nn) = i nn = nn-1 end do end do endif if ( LSAME(cpy_from,'CCS') ) then *--- expand ja indices ... nn = nz do j = ncol,1,-1 do k = 1,ja(j+1)-ja(j) ja(nn) = j nn = nn-1 end do end do endif * *-- if requested, insert diagonal elements even if they are zero... * if ( cpy_diag.eq.'D' ) then nn = MIN( nrow,ncol ) do i = 1,nn iwsp(i) = 0 end do do k = 1,nz if ( ia(k).eq.ja(k) ) iwsp(ia(k)) = 1 end do do i = 1,nn if ( iwsp(i).eq.0 ) then nz = nz + 1 ia(nz) = i ja(nz) = i a(nz) = 0.0d0 endif end do endif *--- COO convertion ... if ( LSAME(cpy_to,'COO') ) return *--- CRS convertion ... if ( LSAME(cpy_to,'CRS') ) call zcmpac( nrow,nz,ia,ja,a, iwsp ) *--- CCS convertion ... if ( LSAME(cpy_to,'CCS') ) call zcmpac( ncol,nz,ja,ia,a, iwsp ) END subroutine zgcoov ( x, y ) c*********************************************************************72 c cc ZGCOOV computes y = A * x, with A stored in COO storage format. c implicit none complex*16 x(*), y(*) * *--- Computes y = A*x. A is passed via a fortran `common statement'. *--- A is assumed here to be under the COOrdinates storage format. * *----------------------------------------------------------------------| *-------------------------------NOTE-----------------------------------* * This is an accessory to Expokit and it is not intended to be * * complete. It is supplied primarily to ensure an unconstrained * * distribution and portability of the package. The matrix-vector * * multiplication routines supplied here fit the non symmetric * * storage and for a symmetric matrix, the entire (not half) matrix * * is required. If the sparsity pattern is known a priori, it is * * recommended to use the most advantageous format and to devise * * the most advantageous matrix-vector multiplication routine. * *----------------------------------------------------------------------* *----------------------------------------------------------------------* integer n, nz, nzmax parameter( nzmax = 50000 ) integer ia(nzmax), ja(nzmax) complex*16 a(nzmax) common /CMAT/ a, ia, ja, nz, n integer i, j complex*16 ZERO parameter( ZERO=(0.0d0,0.0d0) ) do j = 1,n y(j) = ZERO end do do i = 1,nz y(ia(i)) = y(ia(i)) + a(i)*x(ja(i)) end do END subroutine zgcrsv ( x, y ) c*********************************************************************72 c cc ZGCRSV computes y = A * x, with A stored in compressed row storage (CRS). c implicit none complex*16 x(*), y(*) * *--- Computes y = A*x. A is passed via a fortran `common statement'. *--- A is assumed to be under the Compress Row Storage (CRS) format. * integer n, nz, nzmax parameter( nzmax = 50000 ) integer ia(nzmax), ja(nzmax) complex*16 a(nzmax) common /CMAT/ a, ia, ja, nz, n integer i, j complex*16 ZERO parameter( ZERO=(0.0d0,0.0d0) ) do i = 1,n y(i) = ZERO do j = ia(i),ia(i+1)-1 y(i) = y(i) + a(j)*x(ja(j)) end do end do END subroutine zgefa ( a, lda, n, ipvt, info ) c*********************************************************************72 c cc ZGEFA factors a complex general linear system.c c implicit none integer lda,n,ipvt(1),info complex*16 a(lda,1) c c zgefa factors a complex*16 matrix by gaussian elimination. c c zgefa is usually called by zgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for zgeco) = (1 + 9/n)*(time for zgefa) . c c on entry c c a complex*16(lda, n) c the matrix to be factored. 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 and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that zgesl or zgedi will divide by zero c if called. use rcond in zgeco for a reliable c indication of singularity. 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 zaxpy,zscal,izamax c fortran dabs c c internal variables c complex*16 t integer izamax,j,k,kp1,l,nm1 c complex*16 zdum double precision cabs1 double precision dreal,dimag complex*16 zdumr,zdumi dreal(zdumr) = zdumr dimag(zdumi) = (0.0d0,-1.0d0)*zdumi cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum)) c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = izamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (cabs1(a(l,k)) .eq. 0.0d0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -(1.0d0,0.0d0)/a(k,k) call zscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (cabs1(a(n,n)) .eq. 0.0d0) info = n return end subroutine zgemm ( transa, transb, m, n, k, alpha, a, lda, b, ldb, $ beta, c, ldc ) c*********************************************************************72 c cc ZGEMM performs a matrix-matrix operation of the form C=alpha*A*B+beta*C. c implicit none CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC COMPLEX*16 ALPHA, BETA * .. Array Arguments .. COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * ZGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X' or op( X ) = conjg( X' ), * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = conjg( A' ). * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = conjg( B' ). * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - COMPLEX*16 . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - COMPLEX*16 array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. Local Scalars .. LOGICAL CONJA, CONJB, NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB COMPLEX*16 TEMP * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * conjugated or transposed, set CONJA and CONJB as true if A and * B respectively are to be transposed but not conjugated and set * NROWA, NCOLA and NROWB as the number of rows and columns of A * and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) CONJA = LSAME( TRANSA, 'C' ) CONJB = LSAME( TRANSB, 'C' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.CONJA ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.CONJB ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE IF( CONJA )THEN * * Form C := alpha*conjg( A' )*B + beta*C. * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + DCONJG( A( L, I ) )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 150, J = 1, N DO 140, I = 1, M TEMP = ZERO DO 130, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 130 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 140 CONTINUE 150 CONTINUE END IF ELSE IF( NOTA )THEN IF( CONJB )THEN * * Form C := alpha*A*conjg( B' ) + beta*C. * DO 200, J = 1, N IF( BETA.EQ.ZERO )THEN DO 160, I = 1, M C( I, J ) = ZERO 160 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 170, I = 1, M C( I, J ) = BETA*C( I, J ) 170 CONTINUE END IF DO 190, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( B( J, L ) ) DO 180, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 180 CONTINUE END IF 190 CONTINUE 200 CONTINUE ELSE * * Form C := alpha*A*B' + beta*C * DO 250, J = 1, N IF( BETA.EQ.ZERO )THEN DO 210, I = 1, M C( I, J ) = ZERO 210 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 220, I = 1, M C( I, J ) = BETA*C( I, J ) 220 CONTINUE END IF DO 240, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 230, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 230 CONTINUE END IF 240 CONTINUE 250 CONTINUE END IF ELSE IF( CONJA )THEN IF( CONJB )THEN * * Form C := alpha*conjg( A' )*conjg( B' ) + beta*C. * DO 280, J = 1, N DO 270, I = 1, M TEMP = ZERO DO 260, L = 1, K TEMP = TEMP + $ DCONJG( A( L, I ) )*DCONJG( B( J, L ) ) 260 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 270 CONTINUE 280 CONTINUE ELSE * * Form C := alpha*conjg( A' )*B' + beta*C * DO 310, J = 1, N DO 300, I = 1, M TEMP = ZERO DO 290, L = 1, K TEMP = TEMP + DCONJG( A( L, I ) )*B( J, L ) 290 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 300 CONTINUE 310 CONTINUE END IF ELSE IF( CONJB )THEN * * Form C := alpha*A'*conjg( B' ) + beta*C * DO 340, J = 1, N DO 330, I = 1, M TEMP = ZERO DO 320, L = 1, K TEMP = TEMP + A( L, I )*DCONJG( B( J, L ) ) 320 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 330 CONTINUE 340 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 370, J = 1, N DO 360, I = 1, M TEMP = ZERO DO 350, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 350 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 360 CONTINUE 370 CONTINUE END IF END IF * RETURN * * End of ZGEMM . * END subroutine zgemv ( trans, m, n, alpha, a, lda, x, incx, $ beta, y, incy ) c*********************************************************************72 c cc ZGEMV performs a matrix-vector operation of the form y=alpha*A*x+beta*y. c implicit none COMPLEX*16 ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or * * y := alpha*conjg( A' )*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*conjg( A' )*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - COMPLEX*16 array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - COMPLEX*16 . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - COMPLEX*16 array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY LOGICAL NOCONJ * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * NOCONJ = LSAME( TRANS, 'T' ) * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y or y := alpha*conjg( A' )*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 110, J = 1, N TEMP = ZERO IF( NOCONJ )THEN DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE ELSE DO 100, I = 1, M TEMP = TEMP + DCONJG( A( I, J ) )*X( I ) 100 CONTINUE END IF Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 110 CONTINUE ELSE DO 140, J = 1, N TEMP = ZERO IX = KX IF( NOCONJ )THEN DO 120, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 120 CONTINUE ELSE DO 130, I = 1, M TEMP = TEMP + DCONJG( A( I, J ) )*X( IX ) IX = IX + INCX 130 CONTINUE END IF Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 140 CONTINUE END IF END IF * RETURN END subroutine zgesl ( a, lda, n, ipvt, b, job ) c*********************************************************************72 c cc ZGESL solves a factored complex general linear system. c implicit none integer lda,n,ipvt(1),job complex*16 a(lda,1),b(1) c c zgesl solves the complex*16 system c a * x = b or ctrans(a) * x = b c using the factors computed by zgeco or zgefa. c c on entry c c a complex*16(lda, n) c the output from zgeco or zgefa. 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 ipvt integer(n) c the pivot vector from zgeco or zgefa. c c b complex*16(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve ctrans(a)*x = b where c ctrans(a) is the conjugate transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if zgeco has set rcond .gt. 0.0 c or zgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call zgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call zgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue 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 zaxpy,zdotc c fortran dconjg c c internal variables c complex*16 zdotc,t integer k,kb,l,nm1 double precision dreal,dimag complex*16 zdumr,zdumi dreal(zdumr) = zdumr dimag(zdumi) = (0.0d0,-1.0d0)*zdumi c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call zaxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call zaxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve ctrans(a) * x = b c first solve ctrans(u)*y = b c do 60 k = 1, n t = zdotc(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/dconjg(a(k,k)) 60 continue c c now solve ctrans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + zdotc(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine zgesv ( n, m, a,lda, ipiv, b,ldb, iflag ) c*********************************************************************72 c cc ZGESV solves multiple general linear systems of the form A*X=B. c implicit none integer N, M, LDA, LDB, IPIV(N), IFLAG complex*16 A(LDA,N), B(LDB,M) integer j call ZGEFA( A,LDA, N, IPIV, IFLAG ) if ( IFLAG.ne.0 ) stop "Error in ZGESV (LU factorisation)" do j = 1,M call ZGESL( A,LDA, N, IPIV,B(1,j), 0 ) end do end subroutine zgexpv ( n, m, t, v, w, tol, anorm, . wsp,lwsp, iwsp,liwsp, matvec, itrace,iflag ) c*********************************************************************72 c cc ZGEXPV computes w = exp(t*A)*v. c implicit none integer n, m, lwsp, liwsp, itrace, iflag, iwsp(liwsp) double precision t, tol, anorm complex*16 v(n), w(n), wsp(lwsp) external matvec *-----Purpose----------------------------------------------------------| * *--- ZGEXPV computes w = exp(t*A)*v * for a Zomplex (i.e., complex double precision) matrix A * * It does not compute the matrix exponential in isolation but * instead, it computes directly the action of the exponential * operator on the operand vector. This way of doing so allows * for addressing large sparse problems. * * The method used is based on Krylov subspace projection * techniques and the matrix under consideration interacts only * via the external routine `matvec' performing the matrix-vector * product (matrix-free method). * *-----Arguments--------------------------------------------------------| * * n : (input) order of the principal matrix A. * * m : (input) maximum size for the Krylov basis. * * t : (input) time at wich the solution is needed (can be < 0). * * v(n) : (input) given operand vector. * * w(n) : (output) computed approximation of exp(t*A)*v. * * tol : (input/output) the requested accuracy tolerance on w. * If on input tol=0.0d0 or tol is too small (tol.le.eps) * the internal value sqrt(eps) is used, and tol is set to * sqrt(eps) on output (`eps' denotes the machine epsilon). * (`Happy breakdown' is assumed if h(j+1,j) .le. anorm*tol) * * anorm : (input) an approximation of some norm of A. * * wsp(lwsp): (workspace) lwsp .ge. n*(m+1)+n+(m+2)^2+4*(m+2)^2+ideg+1 * +---------+-------+---------------+ * (actually, ideg=6) V H wsp for PADE * * iwsp(liwsp): (workspace) liwsp .ge. m+2 * * matvec : external routine for matrix-vector multiplication. * synopsis: matvec( x, y ) * complex*16 x(*), y(*) * computes: y(1:n) <- A*x(1:n) * where A is the principal matrix. * * itrace : (input) running mode. 0=silent, 1=print step-by-step info * * iflag : (output) exit flag. * <0 - bad input arguments * 0 - no problem * 1 - maximum number of steps reached without convergence * 2 - requested tolerance was too high * *-----Accounts on the computation--------------------------------------| * Upon exit, an interested user may retrieve accounts on the * computations. They are located in the workspace arrays wsp and * iwsp as indicated below: * * location mnemonic description * -----------------------------------------------------------------| * iwsp(1) = nmult, number of matrix-vector multiplications used * iwsp(2) = nexph, number of Hessenberg matrix exponential evaluated * iwsp(3) = nscale, number of repeated squaring involved in Pade * iwsp(4) = nstep, number of integration steps used up to completion * iwsp(5) = nreject, number of rejected step-sizes * iwsp(6) = ibrkflag, set to 1 if `happy breakdown' and 0 otherwise * iwsp(7) = mbrkdwn, if `happy brkdown', basis-size when it occured * -----------------------------------------------------------------| * wsp(1) = step_min, minimum step-size used during integration * wsp(2) = step_max, maximum step-size used during integration * wsp(3) = x_round, maximum among all roundoff errors (lower bound) * wsp(4) = s_round, sum of roundoff errors (lower bound) * wsp(5) = x_error, maximum among all local truncation errors * wsp(6) = s_error, global sum of local truncation errors * wsp(7) = tbrkdwn, if `happy breakdown', time when it occured * wsp(8) = t_now, integration domain successfully covered * wsp(9) = hump, i.e., max||exp(sA)||, s in [0,t] (or [t,0] if t<0) * wsp(10) = ||w||/||v||, scaled norm of the solution w. * -----------------------------------------------------------------| * The `hump' is a measure of the conditioning of the problem. The * matrix exponential is well-conditioned if hump = 1, whereas it is * poorly-conditioned if hump >> 1. However the solution can still be * relatively fairly accurate even when the hump is large (the hump * is an upper bound), especially when the hump and the scaled norm * of w [this is also computed and returned in wsp(10)] are of the * same order of magnitude (further details in reference below). * *----------------------------------------------------------------------| *-----The following parameters may also be adjusted herein-------------| * integer mxstep, mxreject, ideg double precision delta, gamma parameter( mxstep = 500, . mxreject = 0, . ideg = 6, . delta = 1.2d0, . gamma = 0.9d0 ) * mxstep : maximum allowable number of integration steps. * The value 0 means an infinite number of steps. * * mxreject: maximum allowable number of rejections at each step. * The value 0 means an infinite number of rejections. * * ideg : the Pade approximation of type (ideg,ideg) is used as * an approximation to exp(H). The value 0 switches to the * uniform rational Chebyshev approximation of type (14,14) * * delta : local truncation error `safety factor' * * gamma : stepsize `shrinking factor' * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * EXPOKIT: Software Package for Computing Matrix Exponentials. * ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 *----------------------------------------------------------------------| complex*16 ZERO, ONE parameter( ZERO=(0.0d0,0.0d0), ONE=(1.0d0,0.0d0) ) integer i, j, k1, mh, mx, iv, ih, j1v, ns, ifree, lfree, iexph, . ireject,ibrkflag,mbrkdwn, nmult, nreject, nexph, nscale, . nstep double precision sgn, t_out, tbrkdwn, step_min,step_max, err_loc, . s_error, x_error, t_now, t_new, t_step, t_old, . xm, beta, break_tol, p1, p2, p3, eps, rndoff, . vnorm, avnorm, hj1j, hump, SQR1 complex*16 hij intrinsic AINT,ABS,CMPLX,DBLE,INT,LOG10,MAX,MIN,NINT,SIGN,SQRT complex*16 ZDOTC double precision DZNRM2 * *--- check restrictions on input parameters ... * iflag = 0 if ( lwsp.lt.n*(m+2)+5*(m+2)**2+ideg+1 ) iflag = -1 if ( liwsp.lt.m+2 ) iflag = -2 if ( m.ge.n .or. m.le.0 ) iflag = -3 if ( iflag.ne.0 ) stop 'bad sizes (in input of ZGEXPV)' * *--- initialisations ... * k1 = 2 mh = m + 2 iv = 1 ih = iv + n*(m+1) + n ifree = ih + mh*mh lfree = lwsp - ifree + 1 ibrkflag = 0 mbrkdwn = m nmult = 0 nreject = 0 nexph = 0 nscale = 0 t_out = ABS( t ) tbrkdwn = 0.0d0 step_min = t_out step_max = 0.0d0 nstep = 0 s_error = 0.0d0 x_error = 0.0d0 t_now = 0.0d0 t_new = 0.0d0 p1 = 4.0d0/3.0d0 1 p2 = p1 - 1.0d0 p3 = p2 + p2 + p2 eps = ABS( p3-1.0d0 ) if ( eps.eq.0.0d0 ) go to 1 if ( tol.le.eps ) tol = SQRT( eps ) rndoff = eps*anorm break_tol = 1.0d-7 *>>> break_tol = tol *>>> break_tol = anorm*tol sgn = SIGN( 1.0d0,t ) call ZCOPY( n, v,1, w,1 ) beta = DZNRM2( n, w,1 ) vnorm = beta hump = beta * *--- obtain the very first stepsize ... * SQR1 = SQRT( 0.1d0 ) xm = 1.0d0/DBLE( m ) p2 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p2/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 * *--- step-by-step integration ... * 100 if ( t_now.ge.t_out ) goto 500 nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) p1 = 1.0d0/beta do i = 1,n wsp(iv + i-1) = p1*w(i) end do do i = 1,mh*mh wsp(ih+i-1) = ZERO end do * *--- Arnoldi loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) do i = 1,j hij = ZDOTC( n, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) call ZAXPY( n, -hij, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) wsp(ih+(j-1)*mh+i-1) = hij end do hj1j = DZNRM2( n, wsp(j1v),1 ) *--- if `happy breakdown' go straightforward at the end ... if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = CMPLX( hj1j ) call ZDSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DZNRM2( n, wsp(j1v),1 ) * *--- set 1 for the 2-corrected scheme ... * 300 continue wsp(ih+m*mh+m+1) = ONE * *--- loop while ireject>> break_tol = tol *>>> break_tol = anorm*tol sgn = SIGN( 1.0d0,t ) SQR1 = SQRT( 0.1d0 ) call ZCOPY( n, v,1, w,1 ) * *--- step-by-step integration ... * 100 if ( t_now.ge.t_out ) goto 500 nmult = nmult + 1 call matvec( w, wsp(iv) ) call ZAXPY( n, ONE, u,1, wsp(iv),1 ) beta = DZNRM2( n, wsp(iv),1 ) if ( beta.eq.0.0d0 ) goto 500 call ZDSCAL( n, 1.0d0/beta, wsp(iv),1 ) do i = 1,mh*mh wsp(ih+i-1) = ZERO end do if ( nstep.eq.0 ) then *--- obtain the very first stepsize ... xm = 1.0d0/DBLE( m ) p1 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p1/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 endif nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) * *--- Arnoldi loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) do i = 1,j hij = ZDOTC( n, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) call ZAXPY( n, -hij, wsp(iv+(i-1)*n),1, wsp(j1v),1 ) wsp(ih+(j-1)*mh+i-1) = hij end do hj1j = DZNRM2( n, wsp(j1v),1 ) *--- if `happy breakdown' go straightforward at the end ... if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = CMPLX( hj1j ) call ZDSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DZNRM2( n, wsp(j1v),1 ) * *--- set 1's for the 3-extended scheme ... * 300 continue wsp(ih+mh*mbrkdwn) = ONE wsp(ih+(m-1)*mh+m) = ZERO do i = 1,k1-1 wsp(ih+(m+i)*mh+m+i-1) = ONE end do * *--- loop while ireject> 1. However the solution can still be * relatively fairly accurate even when the hump is large (the hump * is an upper bound), especially when the hump and the scaled norm * of w [this is also computed and returned in wsp(10)] are of the * same order of magnitude (further details in reference below). * *----------------------------------------------------------------------| *-----The following parameters may also be adjusted herein-------------| * integer mxstep, mxreject, ideg double precision delta, gamma parameter( mxstep = 500, . mxreject = 0, . ideg = 6, . delta = 1.2d0, . gamma = 0.9d0 ) * mxstep : maximum allowable number of integration steps. * The value 0 means an infinite number of steps. * * mxreject: maximum allowable number of rejections at each step. * The value 0 means an infinite number of rejections. * * ideg : the Pade approximation of type (ideg,ideg) is used as * an approximation to exp(H). The value 0 switches to the * uniform rational Chebyshev approximation of type (14,14) * * delta : local truncation error `safety factor' * * gamma : stepsize `shrinking factor' * *----------------------------------------------------------------------| * Roger B. Sidje (rbs@maths.uq.edu.au) * EXPOKIT: Software Package for Computing Matrix Exponentials. * ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 *----------------------------------------------------------------------| * complex*16 ZERO, ONE parameter( ZERO=(0.0d0,0.0d0), ONE=(1.0d0,0.0d0) ) integer i, j, k1, mh, mx, iv, ih, j1v, ns, ifree, lfree, iexph, . ireject,ibrkflag,mbrkdwn, nmult, nreject, nexph, nscale, . nstep double precision sgn, t_out, tbrkdwn, step_min,step_max, err_loc, . s_error, x_error, t_now, t_new, t_step, t_old, . xm, beta, break_tol, p1, p2, p3, eps, rndoff, . vnorm, avnorm, hj1j, hump, SQR1 complex*16 hjj intrinsic AINT,ABS,CMPLX,DBLE,INT,LOG10,MAX,MIN,NINT,SIGN,SQRT complex*16 ZDOTC double precision DZNRM2 * *--- check restrictions on input parameters ... * iflag = 0 if ( lwsp.lt.n*(m+2)+5*(m+2)**2+ideg+1 ) iflag = -1 if ( liwsp.lt.m+2 ) iflag = -2 if ( m.ge.n .or. m.le.0 ) iflag = -3 if ( iflag.ne.0 ) stop 'bad sizes (in input of DHEXPV)' * *--- initialisations ... * k1 = 2 mh = m + 2 iv = 1 ih = iv + n*(m+1) + n ifree = ih + mh*mh lfree = lwsp - ifree + 1 ibrkflag = 0 mbrkdwn = m nmult = 0 nreject = 0 nexph = 0 nscale = 0 t_out = ABS( t ) tbrkdwn = 0.0d0 step_min = t_out step_max = 0.0d0 nstep = 0 s_error = 0.0d0 x_error = 0.0d0 t_now = 0.0d0 t_new = 0.0d0 p1 = 4.0d0/3.0d0 1 p2 = p1 - 1.0d0 p3 = p2 + p2 + p2 eps = ABS( p3-1.0d0 ) if ( eps.eq.0.0d0 ) go to 1 if ( tol.le.eps ) tol = SQRT( eps ) rndoff = eps*anorm break_tol = 1.0d-7 *>>> break_tol = tol *>>> break_tol = anorm*tol sgn = SIGN( 1.0d0,t ) call ZCOPY( n, v,1, w,1 ) beta = DZNRM2( n, w,1 ) vnorm = beta hump = beta * *--- obtain the very first stepsize ... * SQR1 = SQRT( 0.1d0 ) xm = 1.0d0/DBLE( m ) p2 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p2/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 * *--- step-by-step integration ... * 100 if ( t_now.ge.t_out ) goto 500 nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) beta = DZNRM2( n, w,1 ) p1 = 1.0d0/beta do i = 1,n wsp(iv + i-1) = p1*w(i) end do do i = 1,mh*mh wsp(ih+i-1) = ZERO end do * *--- Lanczos loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) if ( j.gt.1 ) . call ZAXPY(n,-wsp(ih+(j-1)*mh+j-2),wsp(j1v-2*n),1,wsp(j1v),1) hjj = ZDOTC( n, wsp(j1v-n),1, wsp(j1v),1 ) call ZAXPY( n, -hjj, wsp(j1v-n),1, wsp(j1v),1 ) hj1j = DZNRM2( n, wsp(j1v),1 ) wsp(ih+(j-1)*(mh+1)) = hjj *--- if `happy breakdown' go straightforward at the end ... if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = CMPLX( hj1j ) wsp(ih+j*mh+j-1) = CMPLX( hj1j ) call ZDSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DZNRM2( n, wsp(j1v),1 ) * *--- set 1 for the 2-corrected scheme ... * 300 continue wsp(ih+m*mh+m-1) = ZERO wsp(ih+m*mh+m+1) = ONE * *--- loop while ireject>> break_tol = tol *>>> break_tol = anorm*tol sgn = SIGN( 1.0d0,t ) SQR1 = SQRT( 0.1d0 ) call ZCOPY( n, v,1, w,1 ) * *--- step-by-step integration ... * 100 if ( t_now.ge.t_out ) goto 500 nmult = nmult + 1 call matvec( w, wsp(iv) ) call ZAXPY( n, ONE, u,1, wsp(iv),1 ) beta = DZNRM2( n, wsp(iv),1 ) if ( beta.eq.0.0d0 ) goto 500 call ZDSCAL( n, 1.0d0/beta, wsp(iv),1 ) do i = 1,mh*mh wsp(ih+i-1) = ZERO end do if ( nstep.eq.0 ) then *--- obtain the very first stepsize ... xm = 1.0d0/DBLE( m ) p1 = tol*(((m+1)/2.72D0)**(m+1))*SQRT(2.0D0*3.14D0*(m+1)) t_new = (1.0d0/anorm)*(p1/(4.0d0*beta*anorm))**xm p1 = 10.0d0**(NINT( LOG10( t_new )-SQR1 )-1) t_new = AINT( t_new/p1 + 0.55d0 ) * p1 endif nstep = nstep + 1 t_step = MIN( t_out-t_now, t_new ) * *--- Lanczos loop ... * j1v = iv + n do 200 j = 1,m nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) if ( j.gt.1 ) . call ZAXPY(n,-wsp(ih+(j-1)*mh+j-2),wsp(j1v-2*n),1,wsp(j1v),1) hjj = ZDOTC( n, wsp(j1v-n),1, wsp(j1v),1 ) call ZAXPY( n, -hjj, wsp(j1v-n),1, wsp(j1v),1 ) hj1j = DZNRM2( n, wsp(j1v),1 ) wsp(ih+(j-1)*(mh+1)) = hjj if ( hj1j.le.break_tol ) then print*,'happy breakdown: mbrkdwn =',j,' h =',hj1j k1 = 0 ibrkflag = 1 mbrkdwn = j tbrkdwn = t_now t_step = t_out-t_now goto 300 endif wsp(ih+(j-1)*mh+j) = CMPLX( hj1j ) wsp(ih+j*mh+j-1) = CMPLX( hj1j ) call ZDSCAL( n, 1.0d0/hj1j, wsp(j1v),1 ) j1v = j1v + n 200 continue nmult = nmult + 1 call matvec( wsp(j1v-n), wsp(j1v) ) avnorm = DZNRM2( n, wsp(j1v),1 ) * *--- set 1's for the 3-extended scheme ... * 300 continue wsp(ih+mh*mbrkdwn) = ONE wsp(ih+m*mh+m-1) = ZERO wsp(ih+(m-1)*mh+m) = ZERO do i = 1,k1-1 wsp(ih+(m+i)*mh+m+i-1) = ONE end do * *--- loop while ireject