c----------------------------------------------------------------------- c Coded by Edmond Chow, chow@cs.umn.edu c----------------------------------------------------------------------- program info1 c---------------------------------------------------------------------- c usage info1.ex < HB_file c c where info1 is the executable generated by makefile, HB_file is a c file containing a matrix stored in Harwell-Boeing matrices. c Info1 will then dump the information into the standard output. c c To use with larger matrices, increase nmax and nzmax. c---------------------------------------------------------------------- implicit none integer nmax, nzmax parameter (nmax = 30000, nzmax = 800000) integer ia(nmax+1),ia1(nmax+1),ja(nzmax),ja1(nzmax) real*8 a(nzmax),a1(nzmax),rhs(1) character title*72, type*3, key*8, guesol*2 logical valued c integer job, iin, nrow,ncol,nnz,ierr, nrhs, iout c-------------- data iin /5/, iout/6/ c-------------- job = 2 nrhs = 0 call readmt (nmax,nzmax,job,iin,a,ja,ia, rhs, nrhs, * guesol,nrow,ncol,nnz,title,key,type,ierr) c---- if not readable return if (ierr .ne. 0) then write (iout,100) ierr 100 format(' **ERROR: Unable to read matrix',/, * ' Message returned fom readmt was ierr =',i3) stop endif valued = (job .ge. 2) c------- call dinfo1(ncol,iout,a,ja,ia,valued,title,key,type,a1,ja1,ia1) c--------------------end------------------------------------------------ c----------------------------------------------------------------------- end program chkio c------------------------------------------------------------------c c test suite for Part I : I/O routines. c c tests the following : gen5pt.f, prtmt, readmt, amd pltmt. c c 1) generates a 100 x 100 5pt matrix, c c 2) prints it with a given format in file 'first.mat' c c 3) reads the matrix from 'first.mat' using readmat c c 4) prints it again in file 'second.mat' in a different format c c 5) makes 4 pic files to show the different options of pltmt. c c these are in job0.pic, job01.pic, job10.pic, job11.pic c c coded by Y. Saad, RIACS, 08/31/1989. c c------------------------------------------------------------------c parameter (nxmax = 20, nmx = nxmax*nxmax) implicit real*8 (a-h,o-z) integer ia(nmx),ja(7*nmx),iau(nmx) real*8 a(7*nmx),rhs(3*nmx),al(6) character title*72, key*8, type*3, guesol*2 c----- open statements ---------------- open (unit=7,file='first.mat') open (unit=8,file='second.mat') open (unit=20,file='job00.pic') open (unit=21,file='job01.pic') open (unit=22,file='job10.pic') open (unit=23,file='job11.pic') c c---- dimension of grid c nx = 10 ny = 10 nz = 1 al(1) = 1.0D0 al(2) = 0.0D0 al(3) = 2.3D1 al(4) = 0.4D0 al(5) = 0.0D0 al(6) = 8.2D-2 c c---- generate grid problem. c call gen57pt (nx,ny,nz,al,0,n,a,ja,ia,iau,rhs) c c---- create the Harwell-Boeing matrix. Start by defining title, c and type. them define format and print it. c write (title,9) nx, ny 9 format('Five-point matrix on a square region', * ' using a ',I2,' by ',I2,' grid *SPARSKIT*') key = 'Fivept10' type= 'RSA' ifmt = 5 job = 3 guesol = 'GX' c c define a right hand side of ones, an initial guess of two's c and an exact solution of three's. c do 2 k=1, 3*n rhs(k) = real( 1 + (k-1)/n ) 2 continue c call prtmt (n,n,a,ja,ia,rhs,guesol,title,key,type, 1 ifmt,job,7) c---- read it again in same matrix a, ja, ia nmax = nmx nzmax = 7*nmx do 3 k=1, 3*n rhs(k) = 0.0 3 continue job = 3 c rewind 7 c nrhs = 3*n c call readmt (nmax,nzmax,job,7,a,ja,ia,rhs,nrhs,guesol, 1 nrow,ncol,nnz,title,key,type,ierr) print *, ' ierr = ', ierr, ' nrhs ' , nrhs c c matrix read. print it again in a different format c ifmt = 102 ncol = nrow job = 3 c call prtmt (nrow,ncol,a,ja,ia,rhs,guesol,title,key,type, 1 ifmt,job,8) c c---- print four pic files c mode = 0 do 10 i=1, 2 do 11 j=1, 2 job = (i-1)*10 +j-1 iout = 20+(i-1)*2+j-1 call pltmt (nrow,ncol,mode,ja,ia,title,key,type,job,iout) 11 continue 10 continue c-------- stop end program convdif c----------------------------------------------------------------------- c this driver will generate a finite element matrix for the c convection-diffusion problem c c - Div ( K(x,y) Grad u ) + C grad u = f c u = 0 on boundary c c (Dirichlet boundary conditions). c----------------------------------------------------------------------- c this code will prompt for desired mesh (from 0 to 9, with one being c a user input one) and will general the matrix in harewell-Boeing format c in the file mat.hb. It also generates two post script files, one c showing the pattern of the matrix (mat.ps) and the other showing the c corresponding mesh (msh.ps). c----------------------------------------------------------------------- c the structure and organization of the fem codes follow very closely c that of the book by Noborou Kikuchi (Finite element methods in c mechanics, Cambridge Univ. press, 1986). c----------------------------------------------------------------------- c coded Y. Saad and S. Ma -- this version dated August 11, 1993 c----------------------------------------------------------------------- implicit none integer maxnx, maxnel parameter (maxnx = 8000, maxnel = 15000) real*8 a(7*maxnx),x(maxnx),y(maxnx),f(3*maxnx) integer ijk(3,maxnel),ichild(12,maxnel),iparnts(2,maxnx), * ia(maxnx),ja(7*maxnx),iwk(maxnx),jwk(maxnx),nodcode(maxnx), * iperm(maxnx) c character matfile*20, title*72, munt*2,key*8, type*3 real size c integer iin,node,nx,nelx,iout,ndeg,na,nmesh,nref,nxmax,nelmax,nb, * ii,nxnew, nelxnew,ierr,job,n,ncol,mode,ptitle,ifmt external xyk, funb, func, fung data iin/7/,node/3/,nx/0/,nelx/0/,iout/8/,ndeg/12/, na/3000/ c-------------------------------------------------------------- c choose starting mesh --- c-------------------------------------------------------------- c files for output c open(unit=10,file='mat.hb') open(unit=11,file='msh.ps') open(unit=12,file='mat.ps') c----------------------------------------------------------------------- print *, ' enter chosen mesh ' read (*,*) nmesh if (nmesh .eq. 0) then print *, 'enter input file for initial mesh ' read(*,'(a20)') matfile open (unit=7,file=matfile) endif c----------------------------------------------------------------------- print *, ' Enter the number of refinements desired ' read (*,*) nref call inmesh (nmesh,iin,nx,nelx,node,x,y,nodcode,ijk,iperm) c c ...REFINE THE GRID c nxmax = maxnx nelmax= maxnel nb = 0 do 10 ii = 1, nref c c estimate the number nx and nelx at next refinement level. c call checkref(nx,nelx,ijk,node,nodcode,nb,nxnew,nelxnew) if (nxnew .gt. nxmax .or. nelxnew .gt. nelmax) then print *, ' Was able to do only ', ii-1 ,' refinements' goto 11 endif c c ...if OK refine all elements c call refall(nx,nelx,ijk,node,ndeg,x,y,ichild,iparnts,nodcode, * nxmax, nelmax, ierr) if (ierr .ne. 0) print *, '** ERROR IN REFALL : ierr =',ierr 10 continue 11 continue c----------------------------------------------------------------------- job = 0 c----------------------------------------------------------------------- c assemble the matrix in CSR format c----------------------------------------------------------------------- call assmbo (nx,nelx,node,ijk,nodcode,x,y, * a,ja,ia,f,iwk,jwk,ierr,xyk,funb,func,fung) n = nx c----------------------Harewell-Boeing matrix--------------------------- c---------------------1---------2---------3---------5---------6 c 12345678901234567890123456789012345678901234567890 title='1Sample matrix from SPARSKIT ' key = 'SPARSKIT' type = 'rua' ifmt = 6 job = 2 call prtmt (n,n,a,ja,ia,f,'NN',title,key,type,ifmt,job,10) c----------------------Plot of mesh------------------------------------- c----------------------------------------------------------------------- size = 6.0 munt = 'in' mode = 0 title ='Finite element mesh ' ptitle = 1 call psgrid (nx,ja,ia,x,y,title,ptitle,size,munt,11) c hsize = 5.6 c vsize = 5.6 c xleft = 0.0 c bot = 0.0 c job = 30 c call texgrd(nx,ja,ia,x,y,munt,size,vsize,hsize, c * xleft,bot,job,title,ptitle,ijk,node,nelx,11) c----------------------Plot of matrix-pattern--------------------------- c----------------------------------------------------------------------- size = 5.5 mode = 0 title = 'Assembled Matrix' ptitle = 1 ncol = 0 iout = 12 call pspltm(n,n,mode,ja,ia,title,ptitle,size,munt,ncol,iwk,12) c xleft = 0.00 c bot = 0.70 c job = 0 c call texplt(nx,nx,mode,ja,ia,munt,size,vsize,hsize,xleft,bot, c * job,title,ptitle,ncol,iwk,12) stop c-----end-of-program-convdif-------------------------------------------- c----------------------------------------------------------------------- end c----------------------------------------------------------------------- c----------------------------------------------------------------------- c contains the functions needed for defining the PDE poroblems. c c first for the scalar 5-point and 7-point PDE c----------------------------------------------------------------------- function afun (x,y,z) real*8 afun, x,y,z afun = -1.0d0 return end function bfun (x,y,z) real*8 bfun, x,y,z bfun = -1.0d0 return end function cfun (x,y,z) real*8 cfun, x,y,z cfun = -1.0d0 return end function dfun (x,y,z) real*8 dfun, x,y,z data gamma /100.0/ c dfun = gamma * exp( x * y ) dfun = 10.d0 return end function efun (x,y,z) real*8 efun, x,y,z data gamma /100.0/ c efun = gamma * exp( (- x) * y ) efun = 0.d0 return end function ffun (x,y,z) real*8 ffun, x,y,z ffun = 0.0 return end function gfun (x,y,z) real*8 gfun, x,y,z gfun = 0.0 return end function hfun(x, y, z) real*8 hfun, x, y, z hfun = 0.0 return end function betfun(side, x, y, z) real*8 betfun, x, y, z character*2 side betfun = 1.0 return end function gamfun(side, x, y, z) real*8 gamfun, x, y, z character*2 side if (side.eq.'x2') then gamfun = 5.0 else if (side.eq.'y1') then gamfun = 2.0 else if (side.eq.'y2') then gamfun = 7.0 else gamfun = 0.0 endif return end c----------------------------------------------------------------------- c functions for the block PDE's c----------------------------------------------------------------------- subroutine afunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue coeff((j-1)*nfree+j) = -1.0d0 2 continue return end subroutine bfunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue coeff((j-1)*nfree+j) = -1.0d0 2 continue return end subroutine cfunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue coeff((j-1)*nfree+j) = -1.0d0 2 continue return end subroutine dfunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue 2 continue return end subroutine efunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue 2 continue return end subroutine ffunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue 2 continue return end subroutine gfunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue 2 continue return end c----------------------------------------------------------------------- c The material property function xyk for the c finite element problem c----------------------------------------------------------------------- subroutine xyk(nel,xyke,x,y,ijk,node) implicit real*8 (a-h,o-z) dimension xyke(2,2), x(*), y(*), ijk(node,*) c c this is the identity matrix. c xyke(1,1) = 1.0d0 xyke(2,2) = 1.0d0 xyke(1,2) = 0.0d0 xyke(2,1) = 0.0d0 return end c----------------------------------------------------------------------- c contains the functions needed for defining the PDE poroblems. c c first for the scalar 5-point and 7-point PDE c----------------------------------------------------------------------- function afun (x,y,z) real*8 afun, x,y,z afun = -1.0 return end function bfun (x,y,z) real*8 bfun, x,y,z bfun = -1.0 return end function cfun (x,y,z) real*8 cfun, x,y,z cfun = -1.0d0 return end function dfun (x,y,z) real*8 dfun, x,y,z dfun = 10.d0 return end function efun (x,y,z) real*8 efun, x,y,z efun = 0.0d0 return end function ffun (x,y,z) real*8 ffun, x,y,z ffun = 0.0 return end function gfun (x,y,z) real*8 gfun, x,y,z gfun = 0.0 return end function hfun(x, y, z) real*8 hfun, x, y, z hfun = 0.0 return end function betfun(side, x, y, z) real*8 betfun, x, y, z character*2 side betfun = 1.0 return end function gamfun(side, x, y, z) real*8 gamfun, x, y, z character*2 side if (side.eq.'x2') then gamfun = 5.0 else if (side.eq.'y1') then gamfun = 2.0 else if (side.eq.'y2') then gamfun = 7.0 else gamfun = 0.0 endif return end c----------------------------------------------------------------------- c functions for the block PDE's c----------------------------------------------------------------------- subroutine afunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue coeff((j-1)*nfree+j) = -1.0d0 2 continue return end subroutine bfunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue coeff((j-1)*nfree+j) = -1.0d0 2 continue return end subroutine cfunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue coeff((j-1)*nfree+j) = -1.0d0 2 continue return end subroutine dfunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue 2 continue return end subroutine efunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue 2 continue return end subroutine ffunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue 2 continue return end subroutine gfunbl (nfree,x,y,z,coeff) real*8 x, y, z, coeff(100) do 2 j=1, nfree do 1 i=1, nfree coeff((j-1)*nfree+i) = 0.0d0 1 continue 2 continue return end c----------------------------------------------------------------------- c The material property function xyk for the c finite element problem c----------------------------------------------------------------------- c subroutine xyk(nel,xyke,x,y,ijk,node) c implicit real*8 (a-h,o-z) c dimension xyke(2,2), x(*), y(*), ijk(node,*) cc cc this is the identity matrix. cc c xyke(1,1) = 1.0d0 c xyke(2,2) = 1.0d0 c xyke(1,2) = 0.0d0 c xyke(2,1) = 0.0d0 cc c return c end c subroutine xyk (xyke, x, y) implicit real*8(a-h,o-z) dimension xyke(2,2) xyke(1,1) = 1. xyke(1,1) = exp(x+y) xyke(1,2) = 0. xyke(2,1) = 0. xyke(2,2) = 1. xyke(2,2) = exp(x+y) return end function funb(x,y) implicit real*8(a-h,o-z) funb = 0. funb = 2.5 funb = 2*x return end function func(x,y) implicit real*8(a-h,o-z) func = 0. func = -1.5 func = -5*y return end function fung(x,y) c Right hand side corresponding to the exact solution of c u = exp(x+y)*x*(1.-x)*y*(1.-y) c (That exact solution is defined in the function exact) implicit real*8(a-h,o-z) c fung = 1. + x c fung = 2*y*(1.-y) + 2*x*(1.-x) c fung = 2*y*(1.-y) + 2*x*(1.-x) +2.5*(1-2*x)*y*(1-y) c 1 - 1.5*(1.-2*y)*x*(1-x) r = exp(x+y) fung = r*r*((x*x+3.*x)*y*(1.-y) + 1 (y*y+3.*y)*x*(1.-x)) 2 + r*(r-2.*x)*(x*x+x-1.)*y*(1.-y) 3 + r*(r+5.*y)*(y*y+y-1.)*x*(1.-x) return end function exact(x,y) C Exact Solution. implicit real*8(a-h,o-z) exact = exp(x+y)*x*(1.-x)*y*(1.-y) return end program markov c----------------------------------------------------------------------- c c program to generate a Markov chain matrix (to test eigenvalue routines c or algorithms for singular systems (in which case use I-A )) c the matrix models simple random walk on a triangular grid. c see additional comments in subroutine. c ----- c just compile this segment and link to the rest of sparskit c (uses subroutine prtmt from MATGEN) c will create a matrix in the HARWELL/BOEING format and put it in c the file markov.mat c c----------------------------------------------------------------------- parameter (nmax=5000, nzmax= 4*nmax) real*8 a(nzmax) integer ja(nzmax), ia(nmax+1) c character title*72,key*8,type*3 open (unit=11,file='markov.mat') c c read - in grid size - will not accept too large grids. c write (6,'(17hEnter grid-size: ,$)') read *, m if (m*(m+1) .gt. 2*nmax ) then print *, ' m too large - unable to produce matrix ' stop endif c c call generator. c call markgen (m, n, a, ja, ia) c----------------------------------------------------------------------- title=' Test matrix from SPARSKIT - markov chain model ' key = 'randwk01' type = 'rua' iout = 11 job = 2 ifmt = 10 call prtmt (n, n, a, ja, ia, x,'NN',title, * key, type,ifmt, job, iout) stop end c subroutine markgen (m, n, a, ja, ia) c----------------------------------------------------------------------- c matrix generator for a markov model of a random walk on a triang. grid c----------------------------------------------------------------------- c this subroutine generates a test matrix that models a random c walk on a triangular grid. This test example was used by c G. W. Stewart ["{SRRIT} - a FORTRAN subroutine to calculate the c dominant invariant subspaces of a real matrix", c Tech. report. TR-514, University of Maryland (1978).] and in a few c papers on eigenvalue problems by Y. Saad [see e.g. LAA, vol. 34, c pp. 269-295 (1980) ]. These matrices provide reasonably easy c test problems for eigenvalue algorithms. The transpose of the c matrix is stochastic and so it is known that one is an exact c eigenvalue. One seeks the eigenvector of the transpose associated c with the eigenvalue unity. The problem is to calculate the c steady state probability distribution of the system, which is c the eigevector associated with the eigenvalue one and scaled in c such a way that the sum all the components is equal to one. c----------------------------------------------------------------------- c parameters c------------ c on entry : c---------- c m = integer. number of points in each direction. c c on return: c---------- c n = integer. The dimension of the matrix. (In fact n is known c to be equal to (m(m+1))/2 ) c a, c ja, c ia = the matrix stored in CSR format. c c----------------------------------------------------------------------- c Notes: 1) the code will actually compute the transpose of the c stochastic matrix that contains the transition probibilities. c 2) It should also be possible to have a matrix generator c with an additional parameter (basically redefining `half' below c to be another parameter and changing the rest accordingly, but c this is not as simple as it sounds). This is not likely to provide c any more interesting matrices. c----------------------------------------------------------------------- real*8 a(*), cst, pd, pu, half integer ja(*), ia(*) c----------------------------------------------------------------------- data half/0.5d0/ c cst = half/real(m-1) c c --- ix counts the grid point (natural ordering used), i.e., c --- the row number of the matrix. c ix = 0 jax = 1 ia(1) = jax c c sweep y coordinates c do 20 i=1,m jmax = m-i+1 c c sweep x coordinates c do 10 j=1,jmax ix = ix + 1 if (j .eq. jmax) goto 2 pd = cst*real(i+j-1) c c north c a(jax) = pd if (i.eq. 1) a(jax) = a(jax)+pd ja(jax) = ix + 1 jax = jax+1 c east a(jax) = pd if (j .eq. 1) a(jax) = a(jax)+pd ja(jax) = ix + jmax jax = jax+1 c south 2 pu = half - cst*real(i+j-3) if ( j .gt. 1) then a(jax) = pu ja(jax) = ix-1 jax = jax+1 endif c west if ( i .gt. 1) then a(jax) = pu ja(jax) = ix - jmax - 1 jax = jax+1 endif ia(ix+1) = jax 10 continue 20 continue n = ix return end :::::::::::::: rexp.f :::::::::::::: program exptest c------------------------------------------------------------------- c c Test program for exponential propagator using Arnoldi approach c This main program is a very simple test using diagonal matrices c (Krylov subspace methods are blind to the structure of the matrix c except for symmetry). This provides a good way of testing the c accuracy of the method as well as the error estimates. c c------------------------------------------------------------------- implicit real*8 (a-h,o-z) parameter (nmax = 400, ih0=60, ndmx=20,nzmax = 7*nmax) real*8 a(nzmax), u(ih0*nmax), w(nmax),w1(nmax),x(nmax),y(nmax) integer ioff(10) data iout/6/, a0/0.0/, b0/1.0/, epsmac/1.d-10/, eps /1.d-10/ c c set dimension of matrix c n = 100 c--------------------------- define matrix ----------------------------- c A is a single diagonal matrix (ndiag = 1 and ioff(1) = 0 ) c----------------------------------------------------------------------- ndiag = 1 ioff(1) = 0 c c-------- entries in the diagonal are uniformly distributed. c h = 1.0d0 / real(n+1) do 1 j=1, n a(j) = real(j+1)* h 1 continue c-------- write (6,'(10hEnter tn: ,$)') read (5,*) tn c write (6,'(36hEpsilon (desired relative accuracy): ,$)') read (5,*) eps c------- write (6,'(36h m (= dimension of Krylov subspace): ,$)') read (5,*) m c------- c define initial conditions: chosen so that solution = (1,1,1,1..1)^T c------- do 2 j=1,n w(j) = dexp(a(j)*tn) 2 w1(j) = w(j) c call expprod (n, m, eps, tn, u, w, x, y, a, ioff, ndiag) c print *, ' final answer ' print *, (w(k),k=1,20) c do 4 k=1,n 4 w1(k) = dexp(-a(k)*tn) * w1(k) print *, ' exact solution ' print *, (w1(k),k=1,20) c c---------- computing actual 2-norm of error ------------------ c t = 0.0d0 do 47 k=1,n 47 t = t+ (w1(k)-w(k))**2 t = dsqrt(t / ddot(n, w,1,w,1) ) c write (6,*) ' final error', t c-------------------------------------------------------------- stop end c------- subroutine oped(n,x,y,diag,ioff,ndiag) c====================================================== c this kernel performs a matrix by vector multiplication c for a diagonally structured matrix stored in diagonal c format c====================================================== implicit real*8 (a-h,o-z) real*8 x(n), y(n), diag(n,ndiag) common nope, nmvec integer j, n, ioff(ndiag) CDIR$ IVDEP do 1 j=1, n y(j) = 0.00 1 continue c do 10 j=1,ndiag io = ioff(j) i1=max0(1,1-io) i2=min0(n,n-io) CDIR$ IVDEP do 9 k=i1,i2 y(k) = y(k)+diag(k,j)*x(k+io) 9 continue 10 continue nmvec = nmvec + 1 nope = nope + 2*ndiag*n return end c program fivept c----------------------------------------------------------------------- c main program for generating 5 point and 7-point matrices in the c Harwell-Boeing format. Creates a file with containing a c harwell-boeing matrix. typical session: c user answer are after the colon c Enter nx, ny, nz : 10 10 1 c Filename for matrix: test.mat c output matrix in data file : test.mat c c nz = 1 will create a 2-D problem c c----------------------------------------------------------------------- integer nmx, nxmax parameter (nxmax = 50, nmx = nxmax*nxmax) c implicit none integer ia(nmx),ja(7*nmx),iau(nmx) real*8 a(7*nmx),rhs(nmx),al(6) character title*72, key*8, type*3, matfile*50, guesol*2 c----------------------------------------------------------------------- integer nx, ny, nz, iout, n, ifmt, job write (6,*) ' ' write(6,'(22hEnter nx, ny, nz : ,$)') read (5,*) nx, ny, nz write(6,'(22hFilename for matrix : ,$)') read(5,'(a50)') matfile open (unit=7,file=matfile) c c boundary condition is partly specified here c c al(1) = 1.0D0 c al(2) = 0.0D0 c al(3) = 2.3D1 c al(4) = 0.4D0 c al(5) = 0.0D0 c al(6) = 8.2D-2 al(1) = 0.0D0 al(2) = 0.0D0 al(3) = 0.0D1 al(4) = 0.0D0 al(5) = 0.0D0 al(6) = 0.0D0 c call gen57pt (nx,ny,nz,al,0,n,a,ja,ia,iau,rhs) iout = 7 c c write out the matrix c guesol='NN' title = * ' 5-POINT TEST MATRIX FROM SPARSKIT ' c '123456789012345678901234567890123456789012345678901234567890 type = 'RUA' key ='SC5POINT' C 12345678 ifmt = 15 job = 2 c upper part only?? c call getu (n, a, ja, ia, a, ja, ia) call prtmt (n,n,a,ja,ia,rhs,guesol,title,key,type, 1 ifmt,job,iout) write (6,*) ' output matrix in data file : ', matfile c stop end program bfivept c----------------------------------------------------------------------- c main program for generating BLOCK 5 point and 7-point matrices in the c Harwell-Boeing format. Creates a file with containing a c harwell-boeing matrix. c c max block size = 5 c max number of grid points = 8000 = ( nx * ny * nz .le. 8000) c matrix dimension = (nx*ny*nz* Block-size**2) .le. 8000 * 25= 200,000 c c typical session: c Enter nx, ny, nz : 10 10 1 c enter block-size : 4 c enter filename for matrix: test.mat c output matrix in data file : test.mat c c nz =1 will create a 2-D problem c----------------------------------------------------------------------- parameter (nxmax = 20, nmx = nxmax*nxmax*nxmax, ntot=nmx*25) integer ia(ntot),ja(ntot),iau(ntot), iao(ntot),jao(ntot) real*8 stencil(7,100), a(ntot), ao(ntot) character title*72,key*8,type*3, matfile*50, guesol*2 c----------------------------------------------------------------------- write (6,*) ' ' write(6,'(22hEnter nx, ny, nz : ,$)') read (5,*) nx, ny, nz write(6,'(22hnfree (Block size) : ,$)') read (5,*) nfree write(6,'(22hFilename for matrix : ,$)') read(5,'(a50)') matfile open (unit=7,file=matfile) c write (6,*) ' output in data file : ', matfile c------------------------------------------------------ na = nfree*nfree c call gen57bl (nx,ny,nz,nfree,na,n,a,ja,ia,iau,stencil) c------------------------------------------------------ print *, ' n=', n, ' nfree ', nfree, ' na =', na call bsrcsr(1,n,nfree, na, a, ja, ia, ao, jao, iao) n = n * nfree ! Apr. 21, 1995 guesol='NN' title = * ' BLOCK 5-POINT TEST MATRIX FROM SPARSKIT ' type = 'RUA' key = 'BLOCK5PT' C 12345678 ifmt = 15 job = 2 iout = 7 call prtmt (n,n,ao,jao,iao,rhs,guesol,title,key,type, 1 ifmt,job,iout) print *, ' output in data file : ', matfile c stop end program rilut c----------------------------------------------------------------------- c test program for ilut preconditioned gmres. c this program generates a sparse matrix using c matgen and then solves a linear system with an c artificial rhs. c----------------------------------------------------------------------- implicit none c integer nmax, nzmax parameter (nmax=5000,nzmax=100000) integer ia(nmax),ja(nzmax),jau(nzmax),ju(nzmax),iw(nmax*3), & iperm(nmax*2),ipar(16), levs(nzmax) real*8 a(nzmax),x(nmax),y(nmax),au(nzmax),vv(nmax,20), * xran(nmax),rhs(nmax),al(nmax),fpar(16) c real t(2), t1, etime c integer nx,ny,nz,n,j,k,ierr,meth,lfil,nwk,im,maxits,iout real*8 tol,permtol,eps,alph,gammax,gammay,alpha external gmres c common /func/ gammax, gammay, alpha c----------------------------------------------------------------------- c pde to be discretized is : c--------------------------- c c -Lap u + gammax exp (xy)delx u + gammay exp (-xy) dely u +alpha u c c where Lap = 2-D laplacean, delx = part. der. wrt x, c dely = part. der. wrt y. c gammax, gammay, and alpha are passed via the commun func. c c----------------------------------------------------------------------- c c data for PDE: c nx = 30 ny = 30 nz = 1 alpha = -50.0 gammax = 10.0 gammay = 10.0 c c data for preconditioner c nwk = nzmax c c data for GMRES c im = 10 eps = 1.0D-07 maxits = 100 iout = 6 permtol = 1.0 ipar(2) = 2 ipar(3) = 2 ipar(4) = 20*nmax ipar(5) = im ipar(6) = maxits fpar(1) = eps fpar(2) = 2.22D-16 c c same initial guess for gmres c c-------------------------------------------------------------- c call gen57 to generate matrix in compressed sparse row format c-------------------------------------------------------------- c c define part of the boundary condition here c al(1) = 0.0 al(2) = 1.0 al(3) = 0.0 al(4) = 0.0 al(5) = 0.0 al(6) = 0.0 call gen57pt(nx,ny,nz,al,0,n,a,ja,ia,ju,rhs) c c zero initial guess to the iterative solvers c do j=1, n xran(j) = 0.d0 enddo print *, 'RILUT: generated a finite difference matrix' print *, ' grid size = ', nx, ' X ', ny, ' X ', nz print *, ' matrix size = ', n c-------------------------------------------------------------- c gnerate right han side = A * (1,1,1,...,1)**T c-------------------------------------------------------------- do k=1,n x(k) = 1.0 enddo call amux(n, x, y, a, ja, ia) c-------------------------------------------------------------- c test all different methods available: c ILU0, MILU0, ILUT and with different values of tol and lfil c ( from cheaper to more expensive preconditioners) c The more accurate the preconditioner the fewer iterations c are required in pgmres, in general. c c-------------------------------------------------------------- do 200 meth = 1, 15 goto (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) meth 1 continue write (iout,*) ' +++++ ILU(0) Preconditioner ++++ ' c t1 = etime(t) call ilu0 (n, a, ja, ia, au, jau, ju, iw, ierr) c t1 = etime(t) - t1 goto 100 2 continue write (iout,*) ' +++++ MILU(0) Preconditioner ++++ ' c t1 = etime(t) call milu0 (n, a, ja, ia, au, jau, ju, iw, ierr) c t1 = etime(t) - t1 goto 100 3 continue write (iout,*) ' +++++ ILUT Preconditioner ++++ ' write (iout,*) ' +++++ tol = 0.0001, lfil=5 ++++ ' tol = 0.0001 lfil = 5 c c t1 = etime(t) call ilut (n,a,ja,ia,lfil,tol,au,jau,ju,nwk,vv,iw,ierr) c t1 = etime(t) - t1 goto 100 4 continue write (iout,*) ' +++++ ILUT Preconditioner ++++ ' write (iout,*) ' +++++ tol = 0.0001, lfil=10 ++++ ' tol = 0.0001 lfil = 10 c c t1 = etime(t) call ilut (n,a,ja,ia,lfil,tol,au,jau,ju,nwk,vv,iw,ierr) c t1 = etime(t) - t1 goto 100 5 continue write (iout,*) ' +++++ ILUT Preconditioner ++++ ' write (iout,*) ' +++++ tol = .0001, lfil=15 ++++ ' tol = 0.0001 lfil = 15 c c t1 = etime(t) call ilut (n,a,ja,ia,lfil,tol,au,jau,ju,nwk,vv,iw,ierr) c t1 = etime(t) - t1 goto 100 6 continue write (iout,*) ' +++++ ILUTP Preconditioner ++++ ' write (iout,*) ' +++++ tol = 0.0001, lfil=5 ++++ ' tol = 0.0001 lfil = 5 c c t1 = etime(t) call ilutp(n,a,ja,ia,lfil,tol,permtol,n,au,jau,ju,nwk, * vv,iw,iperm,ierr) c t1 = etime(t) - t1 goto 100 7 continue write (iout,*) ' +++++ ILUTP Preconditioner ++++ ' write (iout,*) ' +++++ tol = 0.0001, lfil=10 ++++ ' tol = 0.0001 lfil = 10 c c t1 = etime(t) call ilutp(n,a,ja,ia,lfil,tol,permtol,n,au,jau,ju,nwk, * vv,iw,iperm,ierr) c t1 = etime(t) - t1 goto 100 c----------------------------------------------------------------------------- 8 continue write (iout,*) ' +++++ ILUTP Preconditioner ++++ ' write (iout,*) ' +++++ tol = .0001, lfil=15 ++++ ' tol = 0.0001 lfil = 15 c c t1 = etime(t) call ilutp(n,a,ja,ia,lfil,tol,permtol,n,au,jau,ju,nwk, * vv,iw,iperm,ierr) c t1 = etime(t) - t1 goto 100 9 continue write (iout,*) ' +++++ ILUK Preconditioner ++++ ' write (iout,*) ' +++++ lfil=0 ++++ ' lfil = 0 c t1 = etime(t) call iluk(n,a,ja,ia,lfil,au,jau,ju,levs,nwk,vv,iw,ierr) print *, ' nnz for a =', ia(n+1) - ia(1) print *, ' nnz for ilu =', jau(n+1) -jau(1) + n c t1 = etime(t) - t1 goto 100 c 10 continue write (iout,*) ' +++++ ILUK Preconditioner ++++ ' write (iout,*) ' +++++ lfil=1 ++++ ' lfil = 1 c t1 = etime(t) call iluk(n,a,ja,ia,lfil,au,jau,ju,levs,nwk,vv,iw,ierr) print *, ' nnz for a =', ia(n+1) - ia(1) print *, ' nnz for ilu =', jau(n+1) -jau(1) + n c t1 = etime(t) - t1 goto 100 c 11 continue write (iout,*) ' +++++ ILUK Preconditioner ++++ ' write (iout,*) ' +++++ lfil=3 ++++ ' lfil = 3 c t1 = etime(t) call iluk(n,a,ja,ia,lfil,au,jau,ju,levs,nwk,vv,iw,ierr) print *, ' nnz for a =', ia(n+1) - ia(1) print *, ' nnz for ilu =', jau(n+1) -jau(1) + n c t1 = etime(t) - t1 goto 100 c 12 continue write (iout,*) ' +++++ ILUK Preconditioner ++++ ' write (iout,*) ' +++++ lfil=6 ++++ ' lfil = 6 c t1 = etime(t) call iluk(n,a,ja,ia,lfil,au,jau,ju,levs,nwk,vv,iw,ierr) print *, ' nnz for a =', ia(n+1) - ia(1) print *, ' nnz for ilu =', jau(n+1) -jau(1) + n c t1 = etime(t) - t1 goto 100 c 13 continue c----------------------------------------------------------------------- write (iout,*) ' +++++ ILUD Preconditioner ++++ ' write (iout,*) ' +++++ tol=0.075, alpha=0.0 ++++ ' tol = 0.075 alph= 0.0 c t1 = etime(t) call ilud(n,a,ja,ia,alph,tol,au,jau,ju,nwk,vv,iw,ierr) c print *, ' nnz for a =', ia(n+1) - ia(1) print *, ' nnz for ilu =', jau(n+1) -jau(1) + n c t1 = etime(t) - t1 goto 100 c 14 continue write (iout,*) ' +++++ ILUD Preconditioner ++++ ' write (iout,*) ' +++++ tol=0.075, alpha=1.0 ++++ ' tol = 0.075 alph=1.0 c t1 = etime(t) call ilud(n,a,ja,ia,alph,tol,au,jau,ju,nwk,vv,iw,ierr) print *, ' nnz for a =', ia(n+1) - ia(1) print *, ' nnz for ilu =', jau(n+1) -jau(1) + n c t1 = etime(t) - t1 goto 100 15 continue write (iout,*) ' +++++ ILUD Preconditioner ++++ ' write (iout,*) ' +++++ tol=0.01, alpha=1.0 ++++ ' tol = 0.01 c t1 = etime(t) call ilud(n,a,ja,ia,alph,tol,au,jau,ju,nwk,vv,iw,ierr) print *, ' nnz for a =', ia(n+1) - ia(1) print *, ' nnz for ilu =', jau(n+1) -jau(1) + n c t1 = etime(t) - t1 c goto 100 c 100 continue c c check that return was succesful c c print *, ' ILU factorization time ', t1 print *, ' Precon set-up returned with ierr ', ierr if (ierr .ne. 0) goto 200 c-------------------------------------------------------------- c call GMRES c-------------------------------------------------------------- call runrc(n,y,x,ipar,fpar,vv,xran,a,ja,ia,au,jau,ju, + gmres) print *, 'GMRES return status = ', ipar(1) write (iout,*) ' ' 200 continue c c---------------------------------------------------------- c write (iout,*) ' **** SOLUTION **** ' write(iout, 111) (x(k),k=1,n) 111 format (5d15.5) stop end program riters c----------------------------------------------------------------------- c test program for iters -- the basic iterative solvers c c this program generates a sparse matrix using c GEN57PT and then solves a linear system with an c artificial rhs (the solution is a vector of (1,1,...,1)^T). c----------------------------------------------------------------------- c implicit none c implicit real*8 (a-h,o-z) integer nmax, nzmax, maxits,lwk parameter (nmax=5000,nzmax=100000,maxits=60,lwk=nmax*40) integer ia(nmax),ja(nzmax),jau(nzmax),ju(nzmax),iw(nmax*3) integer ipar(16),nx,ny,nz,i,lfil,nwk,nrow,ierr real*8 a(nzmax),sol(nmax),rhs(nmax),au(nzmax),wk(nmax*40) real*8 xran(nmax), fpar(16), al(nmax) real*8 gammax,gammay,alpha,tol external gen57pt,cg,bcg,dbcg,bcgstab,tfqmr,gmres,fgmres,dqgmres external cgnr, fom, runrc, ilut c common /func/ gammax, gammay, alpha c----------------------------------------------------------------------- c pde to be discretized is : c--------------------------- c c -Lap u + gammax exp (xy)delx u + gammay exp (-xy) dely u +alpha u c c where Lap = 2-D laplacean, delx = part. der. wrt x, c dely = part. der. wrt y. c gammax, gammay, and alpha are passed via the commun func. c c----------------------------------------------------------------------- c c data for PDE: c nx = 6 ny = 6 nz = 1 alpha = 0.0 gammax = 0.0 gammay = 0.0 c c set the parameters for the iterative solvers c ipar(2) = 2 ipar(3) = 1 ipar(4) = lwk ipar(5) = 10 ipar(6) = maxits fpar(1) = 1.0D-5 fpar(2) = 1.0D-10 c-------------------------------------------------------------- c call GEN57PT to generate matrix in compressed sparse row format c c al(1:6) are used to store part of the boundary conditions c (see documentation on GEN57PT.) c-------------------------------------------------------------- al(1) = 0.0 al(2) = 0.0 al(3) = 0.0 al(4) = 0.0 al(5) = 0.0 al(6) = 0.0 nrow = nx * ny * nz call gen57pt(nx,ny,nz,al,0,nrow,a,ja,ia,ju,rhs) print *, 'RITERS: generated a finite difference matrix' print *, ' grid size = ', nx, ' X ', ny, ' X ', nz print *, ' matrix size = ', nrow c c set-up the preconditioner ILUT(15, 1E-4) ! new definition of lfil c lfil = 3 tol = 1.0D-4 nwk = nzmax call ilut (nrow,a,ja,ia,lfil,tol,au,jau,ju,nwk, * wk,iw,ierr) ipar(2) = 2 c c generate a linear system with known solution c do i = 1, nrow sol(i) = 1.0D0 xran(i) = 0.d0 end do call amux(nrow, sol, rhs, a, ja, ia) print *, ' ' print *, ' *** CG ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + cg) print *, ' ' print *, ' *** BCG ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + bcg) print *, ' ' print *, ' *** DBCG ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + dbcg) print *, ' ' print *, ' *** CGNR ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + cgnr) print *, ' ' print *, ' *** BCGSTAB ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + bcgstab) print *, ' ' print *, ' *** TFQMR ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + tfqmr) print *, ' ' print *, ' *** FOM ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + fom) print *, ' ' print *, ' *** GMRES ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + gmres) print *, ' ' print *, ' *** FGMRES ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + fgmres) print *, ' ' print *, ' *** DQGMRES ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + dqgmres) stop end c-----end-of-main c----------------------------------------------------------------------- program riter2 c----------------------------------------------------------------------- c test program for iters -- the basic iterative solvers c c this program reads a Harwell/Boeing matrix from standard input c and solves the linear system with an artifical right-hand side c (the solution is a vector of (1,1,...,1)^T) c----------------------------------------------------------------------- c implicit none c implicit real*8 (a-h,o-z) integer nmax, nzmax, maxits,lwk parameter (nmax=5000,nzmax=100000,maxits=60,lwk=nmax*40) integer ia(nmax),ja(nzmax),jau(nzmax),ju(nzmax),iw(nmax*3) integer ipar(16),i,lfil,nwk,nrow,ierr real*8 a(nzmax),sol(nmax),rhs(nmax),au(nzmax),wk(nmax*40) real*8 xran(nmax), fpar(16), tol character guesol*2, title*72, key*8, type*3 external cg,bcg,dbcg,bcgstab,tfqmr,gmres,fgmres,dqgmres external cgnr, fom, runrc, ilut c c set the parameters for the iterative solvers c ipar(2) = 2 ipar(3) = 1 ipar(4) = lwk ipar(5) = 16 ipar(6) = maxits fpar(1) = 1.0D-5 fpar(2) = 1.0D-10 c-------------------------------------------------------------- c read in a matrix from standard input c-------------------------------------------------------------- iounit = 5 job = 2 nrhs = 0 call readmt (nmax,nzmax,job,iounit,a,ja,ia,a,nrhs, * guesol,nrow,ncol,nnz,title,key,type,ierr) print *, 'READ the matrix ', key, type print *, title print * c c set-up the preconditioner ILUT(15, 1E-4) ! new definition of lfil c lfil = 15 tol = 1.0D-4 ! this is too high for ilut for saylr1 tol = 1.0D-7 nwk = nzmax call ilut (nrow,a,ja,ia,lfil,tol,au,jau,ju,nwk, * wk,iw,ierr) ipar(2) = 2 c c generate a linear system with known solution c do i = 1, nrow sol(i) = 1.0D0 xran(i) = 0.D0 end do call amux(nrow, sol, rhs, a, ja, ia) print *, ' ' print *, ' *** CG ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + cg) print *, ' ' print *, ' *** BCG ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + bcg) print *, ' ' print *, ' *** DBCG ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + dbcg) print *, ' ' print *, ' *** CGNR ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + cgnr) print *, ' ' print *, ' *** BCGSTAB ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + bcgstab) print *, ' ' print *, ' *** TFQMR ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + tfqmr) print *, ' ' print *, ' *** FOM ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + fom) print *, ' ' print *, ' *** GMRES ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + gmres) print *, ' ' print *, ' *** FGMRES ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + fgmres) print *, ' ' print *, ' *** DQGMRES ***' call runrc(nrow,rhs,sol,ipar,fpar,wk,xran,a,ja,ia,au,jau,ju, + dqgmres) stop end c-----end-of-main c----------------------------------------------------------------------- program phitest c------------------------------------------------------------------- c Test program for exponential propagator using Arnoldi approach c This main program is a very simple test using diagonal matrices c (Krylov subspace methods are blind to the structure of the matrix c except for symmetry). This provides a good way of testing the c accuracy of the method as well as the error estimates. This test c program tests the phi-function variant instead of the exponential. c c The subroutine phipro which is tested here computes an approximation c to the vector c c w(tn) = w(t0) + tn * phi( - A * tn ) * (r - A w(t0)) c c where phi(z) = (1-exp(z)) / z c c i.e. it solves dw/dt = -A w + r in [t0,t0+ tn] (returns only w(t0+tn)) c c In this program t0=0, tn is input by the user and A is a simple c diagonal matrix. c c------------------------------------------------------------------- implicit real*8 (a-h,o-z) parameter (nmax = 400, ih0=60, ndmx=20,nzmax = 7*nmax) real*8 a(nzmax), u(ih0*nmax), w(nmax),w1(nmax),x(nmax),y(nmax), * r(nmax) integer ioff(10) data iout/6/,a0/0.0/,b0/1.0/,epsmac/1.d-10/,eps/1.d-10/ c c dimendion of matrix c n = 100 c---------------------------define matrix ----------------------------- c A is a single diagonal matrix (ndiag = 1 and ioff(1) = 0 ) c----------------------------------------------------------------------- ndiag = 1 ioff(1) = 0 c c entriesin the diagonal are uniformly distributed. c h = 1.0d0 / real(n+1) do 1 j=1, n a(j) = real(j+1)*h 1 continue c-------- write (6,'(10hEnter tn: ,$)') read (5,*) tn c write (6,'(36hEpsilon (desired relative accuracy): ,$)') read (5,*) eps c------- write (6,'(36h m (= dimension of Krylov subspace): ,$)') read (5,*) m c c define initial conditions to be (1,1,1,1..1)^T c do 2 j=1,n w(j) = 1.0 - 1.0/real(j+10) r(j) = 1.0 w1(j) = w(j) 2 continue c c c HERE IT IS DONE IN 10 SUBSTEPS c nsteps = 10 c tnh = tn/real(nsteps) c MARCHING LOOP c do jj =1, nsteps c call phiprod (n, m, eps, tnh, u, w, r, x, y, a, ioff, ndiag) call phiprod (n, m, eps, tn, u, w, r, x, y, a, ioff, ndiag) c enddo print *, ' final answer ' print *, (w(k),k=1,20) c c exact solution c do 4 k=1,n w1(k)=w1(k)+((1.0 - dexp(-a(k)*tn) )/a(k))*(r(k) + - a(k)*w1(k) ) 4 continue c c exact solution c print *, ' exact answer ' print *, (w1(k),k=1,20) c c computing actual 2-norm of error c t = 0.0d0 do 47 k=1,n t = t+ (w1(k)-w(k))**2 47 continue t = dsqrt(t / ddot(n, w,1,w,1) ) c write (6,*) ' final error', t c----------------------------------------------------------------------- stop end c----------------------------------------------------------------------- subroutine oped(n,x,y,diag,ioff,ndiag) c====================================================== c this kernel performs a matrix by vector multiplication c for a diagonally structured matrix stored in diagonal format c====================================================== implicit real*8 (a-h,o-z) real*8 x(n), y(n), diag(n,ndiag) common nope, nmvec integer j, n, ioff(ndiag) CDIR$ IVDEP do 1 j=1, n y(j) = 0.00 1 continue c do 10 j=1,ndiag io = ioff(j) i1=max0(1,1-io) i2=min0(n,n-io) CDIR$ IVDEP do 9 k=i1,i2 y(k) = y(k)+diag(k,j)*x(k+io) 9 continue 10 continue nmvec = nmvec + 1 nope = nope + 2*ndiag*n return end c program rsobel integer n, ia(1:200), ja(1:1000), ib(1:200), jb(1:1000) integer nrowc, ncolc integer ic(1:200), jc(1:1000), ierr real*8 a(1:1000), b(1:1000), c(1:1000) write (*, '(1x, 9hInput n: ,$)') read *, n call sobel(n,nrowc,ncolc,c,jc,ic,a,ja,ia,b,jb,ib,1000,ierr) print *, 'ierr =', ierr print *, 'Nrow =', nrowc, ' Ncol =', ncolc call dump(1, nrowc, .true., c, jc, ic, 6) end program runall c----------------------------------------------------------------------- c not present: CG and FGMRES c----------------------------------------------------------------------- c program for running all examples -- for book tables. c----------------------------------------------------------------------- parameter (nmax=10000,np1=nmax+1,nt2=2*nmax,nzmax = 100000, * lw=31*nmax,nzumax = 20*nzmax) c----------------------------------------------------------------------- implicit none c integer ia(nt2),ja(nzmax),jlu(nzumax),levs(nzumax),ju(nmax+1) real*8 wk(lw),rhs(nmax),a(nzmax),alu(nzumax) integer jw(nt2), iperm(nt2) real*8 w(nt2),sol(nt2) integer n, mbloc character guesol*2, title*72, key*8, type*3 c----------------------------------------------------------------------- integer iout,lfil,j,k,maxits,i,ncol,numat,mat,iwk,nnz,ipo,ipi,iin, * job,nrhs,im,outf,ipre,its,ierr,ipre1,ipre2,lf,it,itp, * icode,imat,kk,lfinc, niter character*50 filnam integer iters(20) c c up to 10 matrices c real*8 droptol, dropinc, t, permtol, tol, eps, rand, alph, dnrm2 c c----------------------------------------------------------------------- data iin/4/, iout/8/ c----------------------------------------------------------------------- iwk = nzumax c c read in matrix pathname from file. == then read matrix itself c open (2, file ='matfile') read (2,*) numat i=20 c----------------------------------------------------------------------- open (iin,file='inputs') read (iin,*) alph read (iin,*) im, maxits read (iin,*) eps read (iin,*) ipre1, ipre2, niter read (iin,*) tol, dropinc, permtol read (iin,*) lf, lfinc outf = 30 c c BIG MATRIX LOOP c c----------------------------------------------------------------------- c LOOP THROUGH MATRICES c----------------------------------------------------------------------- do 300 mat =1, numat read (2,'(a50)') filnam imat = 20 open (imat,file =filnam) write (7,*) filnam job = 3 nrhs = nmax call readmt (nmax,nzmax,job,imat,a,ja,ia, rhs, nrhs, * guesol,n,ncol,nnz,title,key,type,ierr) c----------------------------------------------------------------------- rewind(imat) write (7,*) ' -- read -- returned with ierr = ', ierr write (7,*) ' matrix ', filnam write (7,*) ' *** n = ', n, ' nnz ', nnz if (ierr .ne. 0) stop ' not able to read *** ' c c sort matrix c call csort(n,a,ja,ia,jlu,.true.) job = 0 call roscal(n,job,1,a,ja,ia,wk,a,ja,ia,ierr) call coscal(n,job,1,a,ja,ia,wk,a,ja,ia,ierr) c----------------------------------------------------------------------- c LOOP THROUGH PRECONDITONERS c----------------------------------------------------------------------- droptol = tol lfil = lf kk = 0 do 500 ipre = ipre1, ipre2 c----------------------------------------------------------------------- c ipre = 1 --> ILUD c ipre = 2 --> ILUT c ipre = 3 --> ILUDP c ipre = 4 --> ILUTP c ipre = 5 --> ILUK c----------------------------------------------------------------------- lfil = lf droptol = tol do 350 itp =1, niter do k=1,n sol(k) = 1.0 enddo c call amux (n, sol, rhs, a, ja, ia) c goto (1,2,3,4) ipre c 1 continue call ilud (n,a,ja,ia,alph,droptol,alu,jlu,ju,iwk, * w,jw,ierr) if (ierr .ne. 0) write (7,*) ' ilud ierr = ', ierr call stb_test(n,rhs,alu,jlu,ju,t) write (7,*) ' ILUD-drpt =', droptol, ' stbtest = ', t write (7,*) ' --- fill-in =', jlu(n+1)-jlu(1) c goto 100 2 continue c call ilutn (n,a,ja,ia,lfil,droptol,alu,jlu,ju,iwk,w,jw, * ierr) if (ierr .ne. 0) write (7,*) ' ILUT ierr = ', ierr call stb_test(n,rhs,alu,jlu,ju,t) write (7,*) ' ILUT-drpt =', droptol, * ' lfil = ', lfil, ' stbtest = ', t write (7,*) ' --- fill-in =', jlu(n+1)-jlu(1) c goto 100 c c pivoting codes c 3 continue mbloc = n call iludp (n,a,ja,ia,alph,droptol,permtol,mbloc,alu, * jlu,ju,iwk,w,jw,iperm,ierr) if (ierr .ne. 0) write (7,*) ' ilud ierr = ', ierr call stb_test(n,rhs,alu,jlu,ju,t) write (7,*) ' ILUD-drpt =', droptol, ' stbtest = ', t write (7,*) ' --- fill-in =', jlu(n+1)-jlu(1) goto 100 4 continue mbloc = n call ilutpn (n,a,ja,ia,lfil,droptol,permtol,mbloc,alu, * jlu,ju,iwk,w,jw,iperm,ierr) if (ierr .ne. 0) write (7,*) ' milutp ierr = ', ierr c call stb_test(n,rhs,alu,jlu,ju,t) write (7,*) ' ILUTP-drpt =', droptol, * ' lfil = ', lfil, ' stbtest = ', t write (7,*) ' --- fill-in =', jlu(n+1)-jlu(1) goto 100 5 continue c call ilukn(n,a,ja,ia,lfil,alu,jlu,ju,levs,iwk,w,jw,ierr) if (ierr .ne. 0) write (7,*) ' ILUT ierr = ', ierr call stb_test(n,rhs,alu,jlu,ju,t) write (7,*) ' ILUK -- lfil = ', lfil,' stbtest = ', t write (7,*) ' --- fill-in =', jlu(n+1)-jlu(1) c goto 100 c ----------- c big loop -- c ----------- 100 continue do 11 k=1,n sol(k) = 1.0 11 continue call amux (n, sol, rhs, a, ja, ia) do 12 j=1, n sol(j) = rand(j) 12 continue its = 0 c----------------------------------------------------------------------- 10 continue c call fgmr (n,im,rhs,sol,it,wk,ipo,ipi,eps,maxits,0,icode) c----------------------------------------------------------------------- if (icode.eq.1) then call amux(n, wk(ipo), wk(ipi), a, ja, ia) its = its+1 goto 10 else if (icode.eq.3 .or. icode.eq.5) then c c output the residuals here ... c if (ipre .eq. 0) then c NOPRE call dcopy(n, wk(ipo),1,wk(ipi),1) else call lusol(n,wk(ipo),wk(ipi),alu,jlu,ju) endif goto 10 endif c c done *** c call amux(n,sol,wk,a,ja,ia) c do i = 1, n wk(n+i) = sol(i) -1.0D0 wk(i) = wk(i) - rhs(i) end do c kk = kk+1 iters(kk) = its c res(kk) = dnrm2(n,wk,1) c err(kk) = dnrm2(n,wk(n+1),1) write (iout, *) '# actual residual norm ', dnrm2(n,wk,1) write (iout, *) '# error norm ', dnrm2(n,wk(1+n),1) c droptol = droptol*dropinc lfil = lfil + lfinc c c pivoting codes c if (ipre .ge. 3) then do k=1, ia(n+1)-1 ja(k) = iperm(ja(k)) enddo endif 350 continue c 500 continue call entline(outf,filnam,iters,kk) close (imat) 300 continue stop c-------------end-of-main-program-ilut_solve_new------------------------ c----------------------------------------------------------------------- end c----------------------------------------------------------------------- function distdot(n,x,ix,y,iy) integer n, ix, iy real*8 distdot, x(*), y(*), ddot external ddot distdot = ddot(n,x,ix,y,iy) return c-----end-of-distdot c----------------------------------------------------------------------- end c----------------------------------------------------------------------- subroutine entline(outf,mat,its,kk) implicit none integer outf, kk,its(kk),k character mat*70 c real*8 err(kk), res(kk) c----------------------------------------------------------------------- write(outf,100) mat(19:29),(its(k),k=1,8) 100 format * (a,' & ', 8(i4,' & '),'\\\\ \\hline') return end c----------------------------------------------------------------------- subroutine stb_test(n,sol,alu,jlu,ju,tmax) implicit none integer n, jlu(*),ju(*) real*8 sol(n),alu(*),tmax c----------------------------------------------------------------------- real*8 max, t, abs integer j c call lusol(n,sol,sol,alu,jlu,ju) c tmax = 0.0 do j=1, n t = abs(sol(j)) tmax = max(t,tmax) enddo return end program zlatev c----------------------------------------------------------------------- c c test suite for zlatev matrices. generates three matrices and c writes them in three different files in Harwell-Boeing format. c zlatev1.mat produced from matrf2 c zlatev2.mat produced from dcn c zlatev3.mat produced from ecn c c----------------------------------------------------------------------- parameter (nmax = 1000, nzmax=20*nmax) implicit real*8 (a-h,o-z) integer ia(nzmax), ja(nzmax), iwk(nmax) real*8 a(nzmax) character title*72, key*3,type*8, guesol*2 c open (unit=7,file='zlatev1.mat') open (unit=8,file='zlatev2.mat') open (unit=9,file='zlatev3.mat') c m = 100 n = m ic = n/2 index = 10 alpha = 5.0 nn = nzmax c c call matrf2 c call matrf2(m,n,ic,index,alpha,nn,nz,a,ia,ja,ierr) job = 1 c do 110 i = 1, nz c print *, ia(i), ja(i), a(i) c 110 continue call coicsr(n, nz, job, a, ja, ia, iwk) c----- title = ' 1st matrix from zlatev examples ' type = 'RUA' key = ' ZLATEV1' iout = 7 guesol='NN' c ifmt = 3 job = 2 c c write result in H-B format. c c Replaces prtmt with smms in order to print matrix in format for c SMMS instead. c call smms (n,1,n,0,a,ja,ia,iout) call prtmt (n,n,a,ja,ia,rhs,guesol,title,type,key, 1 ifmt,job,iout) c-------- second type of matrices dcn matrices --------------- n = 200 nn = nzmax ic = 20 c------------------------------------------------------- c matrix of the type e(c,n) c------------------------------------------------------- call dcn(a,ia,ja,n,ne,ic,nn,ierr) c--------------------------------------------------- call coicsr(n, ne, job, a, ja, ia, iwk) title = ' 2nd matrix from zlatev examples ' iout = iout+1 guesol='NN' type = 'RUA' key = ' ZLATEV2' c ifmt = 3 job = 2 c c write result in second file c c Replaced prtmt with smms in order to print matrix in format for c SMMS instead. c call smms (n,1,n,0,a,ja,ia,iout) call prtmt (n,n,a,ja,ia,rhs,guesol,title,type,key, 1 ifmt,job,iout) c------------------------------------------------------- c matrix of the type e(c,n) c------------------------------------------------------- n = 200 ic = 20 nn = nzmax c c call ecn c call ecn(n,ic,ne,ia,ja,a,nn,ierr) call coicsr(n, ne, job, a, ja, ia, iwk) title = ' 3nd matrix from zlatev examples ' guesol='NN' type = 'RUA' key = ' ZLATEV3' iout = iout+1 c ifmt = 3 job = 2 c c write resulting matrix in third file c c Replaced prtmt with smms in order to print matrix in format for c SMMS instead. c call smms (n,1,n,0,a,ja,ia,iout) call prtmt (n,n,a,ja,ia,rhs,guesol,title,type,key, 1 ifmt,job,iout) stop end