From 8919dacf4923959e74bc3aa61c5c6ed7c785cdcb Mon Sep 17 00:00:00 2001 From: Alfredo Buttari Date: Thu, 6 Oct 2005 15:56:25 +0000 Subject: [PATCH] cleanup --- test/pargen/RUNS/ppde.inp | 6 +- test/pargen/pp2d.f90 | 640 ---------------------------------- test/pargen/pp2ds.f90 | 700 ------------------------------------- test/pargen/ppde90.f90 | 34 +- test/pargen/ppde90log.f90 | 682 ------------------------------------ test/pargen/ppde90s.f90 | 707 -------------------------------------- 6 files changed, 26 insertions(+), 2743 deletions(-) delete mode 100644 test/pargen/pp2d.f90 delete mode 100644 test/pargen/pp2ds.f90 delete mode 100644 test/pargen/ppde90log.f90 delete mode 100644 test/pargen/ppde90s.f90 diff --git a/test/pargen/RUNS/ppde.inp b/test/pargen/RUNS/ppde.inp index 8c11a1c2..97ee115a 100644 --- a/test/pargen/RUNS/ppde.inp +++ b/test/pargen/RUNS/ppde.inp @@ -1,8 +1,8 @@ 7 Number of entries below this -CGS Iterative method BICGSTAB CGS BICG BICGSTABL -4 Preconditioner ILU DIAGSC NONE +BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL +7 Preconditioner ILU DIAGSC NONE 2 Number ov overlapping levels -CSR A Storage format CSR COO JAD +COO A Storage format CSR COO JAD 20 Domain size (acutal sistem is this**3) 1 Stopping criterion 80 MAXIT diff --git a/test/pargen/pp2d.f90 b/test/pargen/pp2d.f90 deleted file mode 100644 index 207797be..00000000 --- a/test/pargen/pp2d.f90 +++ /dev/null @@ -1,640 +0,0 @@ -! -! This sample program shows how to build and solve a sparse linear -! -! The program solves a linear system based on the partial differential -! equation -! -! -! -! the equation generated is: -! b1 d d (u) b2 d d (u) a1 d (u)) a2 d (u))) -! - ------ - ------ + ----- + ------ + a3 u = 0 -! dx dx dy dy dx dy -! -! -! with Dirichlet boundary conditions on the unit cube -! -! 0<=x,y<=1 -! -! The equation is discretized with finite differences and uniform stepsize; -! the resulting discrete equation is -! -! ( u(x,y)(2b1+2b2+a1+a2)+u(x-1,y)(-b1-a1)+u(x,y-1)(-b2-a2)+ -! -u(x+1,y)b1-u(x,y+1)b2)*(1/h**2) -! -! Example taken from: C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 -! -! -! In this sample program the index space of the discretized -! computational domain is first numbered sequentially in a standard way, -! then the corresponding vector is distributed according to an HPF BLOCK -! distribution directive. -! -! Boundary conditions are set in a very simple way, by adding -! equations of the form -! -! u(x,y) = rhs(x,y) -! -Program PP2D - USE F90SPARSE - Implicit none - - interface - !.....user passed subroutine..... - subroutine part_block(glob_index,n,np,pv,nv) - INTEGER, INTENT(IN) :: GLOB_INDEX, N, NP - INTEGER, INTENT(OUT) :: NV - INTEGER, INTENT(OUT) :: PV(*) - end subroutine part_block - end interface - ! input parameters - Character :: CMETHD*10, PREC*10, AFMT*5 - Integer :: IDIM, IRET - - ! Miscellaneous - Integer, Parameter :: IZERO=0, IONE=1 - Character, PARAMETER :: ORDER='R' - INTEGER :: IARGC,CONVERT_DESCR,dim, CHECK_DESCR - REAL(KIND(1.D0)), PARAMETER :: DZERO = 0.D0, ONE = 1.D0 - REAL(KIND(1.D0)) :: MPI_WTIME, T1, T2, TPREC, TSOLVE, T3, T4 - EXTERNAL MPI_WTIME - - ! Sparse Matrix and preconditioner - TYPE(D_SPMAT) :: A, L, U, H - TYPE(D_PREC) :: PRE - ! Descriptor - TYPE(desc_type) :: DESC_A, DESC_A_OUT - ! Dense Matrices - REAL(KIND(1.d0)), POINTER :: B(:), X(:), D(:),LD(:) - INTEGER, pointer :: WORK(:) - ! BLACS parameters - INTEGER :: nprow, npcol, icontxt, iam, np, myprow, mypcol - - ! Solver parameters - INTEGER :: ITER, ITMAX,IERR,ITRACE, METHD,IPREC, ISTOPC,& - & IPARM(20), ML - REAL(KIND(1.D0)) :: ERR, EPS, RPARM(20) - - ! Other variables - INTEGER :: I,INFO - INTEGER :: INTERNAL, M,II - - ! Initialize BLACS - CALL BLACS_PINFO(IAM, NP) - CALL BLACS_GET(IZERO, IZERO, ICONTXT) - - ! Rectangular Grid, P x 1 - - CALL BLACS_GRIDINIT(ICONTXT, ORDER, NP, IONE) - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - ! - ! Get parameters - ! - CALL GET_PARMS(ICONTXT,CMETHD,PREC,AFMT,IDIM,ISTOPC,ITMAX,ITRACE,ML) - - ! - ! Allocate and fill in the coefficient matrix, RHS and initial guess - ! - - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - CALL CREATE_MATRIX(IDIM,A,B,X,DESC_A,PART_BLOCK,ICONTXT,AFMT) - T2 = MPI_WTIME() - T1 - - DIM=SIZE(A%ASPK) - - ALLOCATE(H%ASPK(DIM),H%IA1(DIM),H%IA2(DIM),H%PL(SIZE(A%PL)),& - & H%PL(SIZE(A%PL)),D(SIZE(A%PL)),& - & DESC_A_OUT%MATRIX_DATA(SIZE(DESC_A%MATRIX_DATA)),& - & DESC_A_OUT%HALO_INDEX(SIZE(DESC_A%HALO_INDEX)),& - & DESC_A_OUT%OVRLAP_INDEX(SIZE(DESC_A%OVRLAP_INDEX)),& - & DESC_A_OUT%OVRLAP_ELEM(SIZE(DESC_A%OVRLAP_ELEM)),& - & DESC_A_OUT%LOC_TO_GLOB(SIZE(DESC_A%LOC_TO_GLOB)),& - & DESC_A_OUT%GLOB_TO_LOC(SIZE(DESC_A%GLOB_TO_LOC)), WORK(1024)) - check_descr=15 -! work(5)=9 -!!$ WRITE(0,*)'CALLING VERIFY' -!!$ CALL F90_PSVERIFY(D,A,DESC_A,CHECK_DESCR,CONVERT_DESCR,H,& -!!$ & DESC_A_OUT,WORK) -!!$ WRITE(0,*)'VERIFY DONE',CONVERT_DESCR - - deallocate(work) - - CALL DGAMX2D(ICONTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1) - IF (IAM.EQ.0) Write(6,*) 'Matrix creation Time : ',T2 - - - ! - ! Prepare the preconditioner. - ! - SELECT CASE (PREC) - CASE ('SCHW') - IPREC = 3 - CASE ('ILU') - IPREC = 2 - CASE ('DIAGSC') - IPREC = 1 - CASE ('NONE') - IPREC = 0 - CASE DEFAULT - WRITE(0,*) 'Unknown preconditioner' - CALL BLACS_ABORT(ICONTXT,-1) - END SELECT - pre%prec=iprec - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - CALL PRECONDITIONER(A,PRE,DESC_A,IRET) -!!$ CALL PRECONDITIONER(IPREC,A,L,U,D,DESC_A,IRET) - TPREC = MPI_WTIME()-T1 - - CALL DGAMX2D(icontxt,'A',' ',IONE, IONE,TPREC,IONE,t1,t1,-1,-1,-1) - - IF (IAM.EQ.0) WRITE(6,*) 'Preconditioner Time : ',TPREC - - IF (IRET.NE.0) THEN - WRITE(0,*) 'Error on preconditioner',IRET - CALL BLACS_ABORT(ICONTXT,-1) - STOP - END IF - - ! - ! Iterative method parameters - ! - write(*,*) 'Calling Iterative method', size(b),ml - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - EPS = 1.D-9 - IF (CMETHD.EQ.'BICGSTAB') THEN - CALL F90_BICGSTAB(A,PRE,B,X,EPS,DESC_A,& - & ITMAX,ITER,ERR,IERR,ITRACE) -!!$ ELSE IF (CMETHD.EQ.'BICG') THEN -!!$ CALL F90_BICG(A,PRE,B,X,EPS,DESC_A,& -!!$ & ITMAX,ITER,ERR,IERR,ITRACE) -!!$ ELSE IF (CMETHD.EQ.'CGS') THEN -!!$ CALL F90_CGS(A,PRE,B,X,EPS,DESC_A,& -!!$ & ITMAX,ITER,ERR,IERR,ITRACE) -!!$ ELSE IF (CMETHD.EQ.'BICGSTABL') THEN -!!$ CALL F90_BICGSTABL(A,PRE,B,X,EPS,DESC_A,& -!!$ & ITMAX,ITER,ERR,IERR,ITRACE,ML) - ELSE - write(0,*) 'Unknown method ',cmethd - end IF - - CALL BLACS_BARRIER(ICONTXT,'All') - T2 = MPI_WTIME() - T1 - CALL DGAMX2D(ICONTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1) - - IF (IAM.EQ.0) THEN - WRITE(6,*) 'Time to Solve Matrix : ',T2 - WRITE(6,*) 'Time per iteration : ',T2/ITER - WRITE(6,*) 'Number of iterations : ',ITER - WRITE(6,*) 'Error on exit : ',ERR - WRITE(6,*) 'INFO on exit : ',IERR - END IF - - ! - ! Cleanup storage and exit - ! - CALL F90_PSDSFREE(B,DESC_A) - CALL F90_PSDSFREE(X,DESC_A) -!!$ CALL F90_PSDSFREE(D,DESC_A) - - CALL F90_PSSPFREE(A,DESC_A) -!!$ CALL F90_PSSPFREE(L,DESC_A) -!!$ CALL F90_PSSPFREE(U,DESC_A) - CALL F90_PSDSCFREE(DESC_A,info) - - CALL BLACS_GRIDEXIT(ICONTXT) - CALL BLACS_EXIT(0) - - STOP - -CONTAINS - ! - ! Get iteration parameters from the command line - ! - SUBROUTINE GET_PARMS(ICONTXT,CMETHD,PREC,AFMT,IDIM,ISTOPC,ITMAX,ITRACE,ML) - integer :: icontxt - Character :: CMETHD*10, PREC*10, AFMT*5 - Integer :: IDIM, IRET, ISTOPC,ITMAX,ITRACE,ML - Character*40 :: CHARBUF - INTEGER :: IARGC, NPROW, NPCOL, MYPROW, MYPCOL - EXTERNAL IARGC - INTEGER :: INTBUF(10), IP - - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - IF (MYPROW==0) THEN - READ(*,*) IP - IF (IP.GE.3) THEN - READ(*,*) CMETHD - READ(*,*) PREC - READ(*,*) AFMT - - ! Convert strings in array - DO I = 1, LEN(CMETHD) - INTBUF(I) = IACHAR(CMETHD(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - DO I = 1, LEN(PREC) - INTBUF(I) = IACHAR(PREC(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - DO I = 1, LEN(AFMT) - INTBUF(I) = IACHAR(AFMT(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - READ(*,*) IDIM - IF (IP.GE.4) THEN - READ(*,*) ISTOPC - ELSE - ISTOPC=1 - ENDIF - IF (IP.GE.5) THEN - READ(*,*) ITMAX - ELSE - ITMAX=500 - ENDIF - IF (IP.GE.6) THEN - READ(*,*) ITRACE - ELSE - ITRACE=-1 - ENDIF - IF (IP.GE.7) THEN - READ(*,*) ML - ELSE - ML=1 - ENDIF - ! Broadcast parameters to all processors - - INTBUF(1) = IDIM - INTBUF(2) = ISTOPC - INTBUF(3) = ITMAX - INTBUF(4) = ITRACE - INTBUF(5) = ML - CALL IGEBS2D(ICONTXT,'ALL',' ',5,1,INTBUF,5) - - WRITE(6,*)'Solving matrix: ELL1' - WRITE(6,*)'on grid',IDIM,'x',IDIM,'x',IDIM - WRITE(6,*)' with BLOCK data distribution, NP=',Np,& - & ' Preconditioner=',PREC,& - & ' Iterative methd=',CMETHD - ELSE - ! Wrong number of parameter, print an error message and exit - CALL PR_USAGE(0) - CALL BLACS_ABORT(ICONTXT,-1) - STOP 1 - ENDIF - ELSE - ! Receive Parameters - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 10 - CMETHD(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 10 - PREC(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 5 - AFMT(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',5,1,INTBUF,5,0,0) - IDIM = INTBUF(1) - ISTOPC = INTBUF(2) - ITMAX = INTBUF(3) - ITRACE = INTBUF(4) - ML = INTBUF(5) - END IF - RETURN - - END SUBROUTINE GET_PARMS - ! - ! Print an error message - ! - SUBROUTINE PR_USAGE(IOUT) - INTEGER :: IOUT - WRITE(IOUT,*)'Incorrect parameter(s) found' - WRITE(IOUT,*)' Usage: pde90 methd prec dim & - &[istop itmax itrace]' - WRITE(IOUT,*)' Where:' - WRITE(IOUT,*)' methd: CGSTAB TFQMR CGS' - WRITE(IOUT,*)' prec : ILU DIAGSC NONE' - WRITE(IOUT,*)' dim number of points along each axis' - WRITE(IOUT,*)' the size of the resulting linear ' - WRITE(IOUT,*)' system is dim**3' - WRITE(IOUT,*)' istop Stopping criterion 1, 2 or 3 [1] ' - WRITE(IOUT,*)' itmax Maximum number of iterations [500] ' - WRITE(IOUT,*)' itrace 0 (no tracing, default) or ' - WRITE(IOUT,*)' >= 0 do tracing every ITRACE' - WRITE(IOUT,*)' iterations ' - END SUBROUTINE PR_USAGE - -! -! Subroutine to allocate and fill in the coefficient matrix and -! the RHS. -! - SUBROUTINE CREATE_MATRIX(IDIM,A,B,T,DESC_A,PARTS,ICONTXT,AFMT) - ! - ! Discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) a1 d(u) a2 d(u) - ! - ------ - ------ - ----- - ------ + a4 u - ! dxdx dydy dx dy - ! - ! = 0 - ! - ! boundary condition: Dirichlet - ! 0< x,y<1 - ! - ! u(x,y)(2b1+2b2+a1+a2)+u(x-1,y)(-b1-a1)+u(x,y-1)(-b2-a2)+ - ! -u(x+1,y,z)b1-u(x,y+1,z)b2 - - USE TYPESP - USE TYPEDESC - USE F90TOOLS - Implicit None - INTEGER :: IDIM -!!$ external parts - integer, parameter :: nbmax=10 - Real(Kind(1.D0)),Pointer :: B(:),T(:) - Type (desc_type) :: DESC_A - Integer :: ICONTXT - INTERFACE - ! .....user passed subroutine..... - SUBROUTINE PARTS(GLOBAL_INDX,N,NP,PV,NV) - INTEGER, INTENT(IN) :: GLOBAL_INDX, N, NP - INTEGER, INTENT(OUT) :: NV, PV(*) - END SUBROUTINE PARTS - END INTERFACE ! Local variables - Type(D_SPMAT) :: A - Real(Kind(1.d0)) :: ZT(NBMAX),GLOB_X,GLOB_Y,GLOB_Z - Integer :: M,N,NNZ,GLOB_ROW,J - Type (D_SPMAT) :: ROW_MAT - Integer :: X,Y,Z,COUNTER,IA,I,INDX_OWNER - INTEGER :: NPROW,NPCOL,MYPROW,MYPCOL - Integer :: ELEMENT - INTEGER :: INFO, NV, INV - INTEGER, ALLOCATABLE :: PRV(:) - INTEGER, pointer :: ierrv(:) - Real(Kind(1.d0)), pointer :: DWORK(:) - INTEGER,POINTER :: IWORK(:) - character :: afmt*5 - ! deltah dimension of each grid cell - ! deltat discretization time - Real(Kind(1.D0)) :: DELTAH - Real(Kind(1.d0)),Parameter :: RHS=0.d0,ONE=1.d0,ZERO=0.d0 - Real(Kind(1.d0)) :: MPI_WTIME, T1, T2, T3, TINS - Real(Kind(1.d0)) :: a1, a2, a3, a4, b1, b2, b3 - external mpi_wtime,a1, a2, a3, a4, b1, b2, b3 - integer :: nb, ir1, ir2, ipr - logical :: own - ! common area - - - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - DELTAH = 1.D0/(IDIM-1) - - ! Initialize array descriptor and sparse matrix storage. Provide an - ! estimate of the number of non zeroes - CALL SETERR(2) - allocate(ierrv(6)) - - ierrv(:) = 0 - M = IDIM*IDIM - N = M - NNZ = ((N*6)/(NPROW*NPCOL)) - write(*,*) 'Size: n ',n - Call F90_PSDSCALL(N,N,PARTS,ICONTXT,IERRV,DESC_A) - write(*,*) 'Allocating A : nnz',nnz - Call F90_PSSPALL(A,IERRV,DESC_A,NNZ=NNZ) - ! Define RHS from boundary conditions; also build initial guess - write(*,*) 'Allocating B' - Call F90_PSDSALL(N,B,IERRV,DESC_A) - write(*,*) 'Allocating T' - Call F90_PSDSALL(N,T,IERRV,DESC_A) - - ! We build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. Might be extended to generate - ! a bunch of rows per call. - ! - ROW_MAT%DESCRA(1:1) = 'G' - ROW_MAT%FIDA = 'CSR' - write(*,*) 'Allocating ROW_MAT',20*nbmax - ALLOCATE(ROW_MAT%ASPK(20*nbmax),ROW_MAT%IA1(20*nbmax),& - &ROW_MAT%IA2(20*nbmax),PRV(NPROW),stat=info) - if (info.ne.0 ) then - write(*,*) 'Memory allocation error' - call blacs_abort(icontxt,-1) - endif - - TINS = 0.D0 - CALL BLACS_BARRIER(ICONTXT,'ALL') - T1 = MPI_WTIME() - - ! Loop over rows belonging to current process in a BLOCK - ! distribution. - - Z=0 - ROW_MAT%IA2(1)=1 - DO GLOB_ROW = 1, N - CALL PARTS(GLOB_ROW,N,NPROW,PRV,NV) - DO INV = 1, NV - INDX_OWNER = PRV(INV) - IF (INDX_OWNER == MYPROW) THEN - ! Local matrix pointer - ELEMENT=1 - ! Compute gridpoint Coordinates - IF (MOD(GLOB_ROW,(IDIM)).EQ.0) THEN - X = GLOB_ROW/(IDIM) - ELSE - X = GLOB_ROW/(IDIM)+1 - ENDIF - Y = GLOB_ROW-(X-1)*IDIM - ! GLOB_X, GLOB_Y, GLOB_X coordinates - GLOB_X=X*DELTAH - GLOB_Y=Y*DELTAH - GLOB_Z=Z*DELTAH - - - ! Check on boundary points - IF (X.EQ.1) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y) - ELEMENT=ELEMENT+1 - ELSE IF (Y.EQ.1) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y) - ELEMENT=ELEMENT+1 - ELSE IF (X.EQ.IDIM) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y) - ELEMENT=ELEMENT+1 - ELSE IF (Y.EQ.IDIM) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y) - ELEMENT=ELEMENT+1 - ELSE - ! Internal point: build discretization - ! - ! Term depending on (x-1,y) - ! - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z)& - & -A1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-2)*IDIM+(Y) - ELEMENT=ELEMENT+1 - ! Term depending on (x,y-1,z) - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z)& - & -A2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y-1) - ELEMENT=ELEMENT+1 -!!$ ! Term depending on (x,y,z-1) -!!$ ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z)& -!!$ & -A3(GLOB_X,GLOB_Y,GLOB_Z) -!!$ ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& -!!$ & DELTAH) -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z-1) -!!$ ELEMENT=ELEMENT+1 - ! Term depending on (x,y,z) - ROW_MAT%ASPK(ELEMENT)=2*B1(GLOB_X,GLOB_Y,GLOB_Z)& - & +2*B2(GLOB_X,GLOB_Y,GLOB_Z)& - & +A1(GLOB_X,GLOB_Y,GLOB_Z)& - & +A2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y) - ELEMENT=ELEMENT+1 -!!$ ! Term depending on (x,y,z+1) -!!$ ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z) -!!$ ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& -!!$ & DELTAH) -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z+1) -!!$ ELEMENT=ELEMENT+1 -!!$ ! Term depending on (x,y+1,z) - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y+1) - ELEMENT=ELEMENT+1 - ! Term depending on (x+1,y,z) - ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X)*IDIM+(Y) - ELEMENT=ELEMENT+1 - ENDIF - ROW_MAT%M=1 - ROW_MAT%K=N - ROW_MAT%IA2(2)=ELEMENT - ! IA== GLOBAL ROW INDEX - IA=GLOB_ROW -!!$ IA=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ write(0,*) 'Inserting row ',ia,' On proc',myprow - T3 = MPI_WTIME() - CALL F90_PSSPINS(A,IA,1,ROW_MAT,IERRV,DESC_A) - if (ierrv(1).ne.0) then - write(0,*) 'On row ',ia,' IERRV:',ierrv(:) - endif - TINS = TINS + (MPI_WTIME()-T3) - ! Build RHS - IF (X==1) THEN - GLOB_Y=(Y-IDIM/2)*DELTAH - GLOB_Z=(Z-IDIM/2)*DELTAH - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2) - ELSE IF ((Y==1).OR.(Y==IDIM).OR.(Z==1).OR.(Z==IDIM)) THEN - GLOB_X=3*(X-1)*DELTAH - GLOB_Y=(Y-IDIM/2)*DELTAH - GLOB_Z=(Z-IDIM/2)*DELTAH - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X) - ELSE - ZT(1) = 0.D0 - ENDIF - CALL F90_PSDSINS(1,B,IA,ZT(1:1),IERRV,DESC_A) - ZT(1)=0.D0 - CALL F90_PSDSINS(1,T,IA,ZT(1:1),IERRV,DESC_A) - END IF - END DO - END DO - - CALL BLACS_BARRIER(ICONTXT,'ALL') - T2 = MPI_WTIME() - - WRITE(*,*) ' pspins time',TINS - WRITE(*,*) ' Insert time',(T2-T1) - - DEALLOCATE(ROW_MAT%ASPK,ROW_MAT%IA1,ROW_MAT%IA2) - - write(*,*) 'Calling SPASB' - CALL BLACS_BARRIER(ICONTXT,'ALL') - T1 = MPI_WTIME() - - CALL F90_PSSPASB(A,IERRV,DESC_A,AFMT=AFMT) - - CALL BLACS_BARRIER(ICONTXT,'ALL') - T2 = MPI_WTIME() - - WRITE(0,*) ' Assembly time',(T2-T1),' ',a%fida(1:4) - - CALL F90_PSDSASB(B,IERRV,DESC_A) - CALL F90_PSDSASB(T,IERRV,DESC_A) - IF (MYPROW.EQ.0) THEN - WRITE(0,*) ' End CREATE_MATRIX' - ENDIF - RETURN - - END SUBROUTINE CREATE_MATRIX -END PROGRAM PP2D -! -! Functions parametrizing the differential equation -! -FUNCTION A1(X,Y,Z) - REAL(KIND(1.D0)) :: A1 - REAL(KIND(1.D0)) :: X,Y,Z - A1=1.D0 -END FUNCTION A1 -FUNCTION A2(X,Y,Z) - REAL(KIND(1.D0)) :: A2 - REAL(KIND(1.D0)) :: X,Y,Z - A2=2.D1*Y -END FUNCTION A2 -FUNCTION A3(X,Y,Z) - REAL(KIND(1.D0)) :: A3 - REAL(KIND(1.D0)) :: X,Y,Z - A3=1.D0 -END FUNCTION A3 -FUNCTION A4(X,Y,Z) - REAL(KIND(1.D0)) :: A4 - REAL(KIND(1.D0)) :: X,Y,Z - A4=1.D0 -END FUNCTION A4 -FUNCTION B1(X,Y,Z) - REAL(KIND(1.D0)) :: B1 - REAL(KIND(1.D0)) :: X,Y,Z - B1=1.D0 -END FUNCTION B1 -FUNCTION B2(X,Y,Z) - REAL(KIND(1.D0)) :: B2 - REAL(KIND(1.D0)) :: X,Y,Z - B2=1.D0 -END FUNCTION B2 -FUNCTION B3(X,Y,Z) - REAL(KIND(1.D0)) :: B3 - REAL(KIND(1.D0)) :: X,Y,Z - B3=1.D0 -END FUNCTION B3 - - diff --git a/test/pargen/pp2ds.f90 b/test/pargen/pp2ds.f90 deleted file mode 100644 index 5ee6610b..00000000 --- a/test/pargen/pp2ds.f90 +++ /dev/null @@ -1,700 +0,0 @@ -! -! This sample program shows how to build and solve a sparse linear -! -! The program solves a linear system based on the partial differential -! equation -! -! -! -! the equation generated is: -! b1 d d (u) b2 d d (u) a1 d (u)) a2 d (u))) -! - ------ - ------ + ----- + ------ + a3 u = 0 -! dx dx dy dy dx dy -! -! -! with Dirichlet boundary conditions on the unit cube -! -! 0<=x,y,z<=1 -! -! The equation is discretized with finite differences and uniform stepsize; -! the resulting discrete equation is -! -! ( u(x,y,z)(2b1+2b2+a1+a2)+u(x-1,y)(-b1-a1)+u(x,y-1)(-b2-a2)+ -! -u(x+1,y)b1-u(x,y+1)b2)*(1/h**2) -! -! Example taken from: C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 -! -! -! In this sample program the index space of the discretized -! computational domain is first numbered sequentially in a standard way, -! then the corresponding vector is distributed according to an HPF BLOCK -! distribution directive. -! -! Boundary conditions are set in a very simple way, by adding -! equations of the form -! -! u(x,y) = rhs(x,y) -! -Program PDE90 - USE F90SPARSE - Implicit none - - interface - !.....user passed subroutine..... - subroutine part_block(glob_index,n,np,pv,nv) - INTEGER, INTENT(IN) :: GLOB_INDEX, N, NP - INTEGER, INTENT(OUT) :: NV - INTEGER, INTENT(OUT) :: PV(*) - end subroutine part_block - end interface - ! input parameters - Character :: CMETHD*10, PREC*10, AFMT*5 - Integer :: IDIM, IRET - - ! Miscellaneous - Integer, Parameter :: IZERO=0, IONE=1 - Character, PARAMETER :: ORDER='R' - INTEGER :: IARGC,CONVERT_DESCR,dim, CHECK_DESCR - REAL(KIND(1.D0)), PARAMETER :: DZERO = 0.D0, ONE = 1.D0 - REAL(KIND(1.D0)) :: MPI_WTIME, T1, T2, TPREC, TSOLVE, T3, T4 - EXTERNAL MPI_WTIME - - ! Sparse Matrix and preconditioner - TYPE(D_SPMAT) :: A, L, U, H - TYPE(D_PREC) :: PRE - ! Descriptor - TYPE(desc_type) :: DESC_A, DESC_A_OUT - ! Dense Matrices - REAL(KIND(1.d0)), POINTER :: B(:), X(:), D(:),LD(:) - INTEGER, pointer :: WORK(:) - ! BLACS parameters - INTEGER :: nprow, npcol, icontxt, iam, np, myprow, mypcol - - ! Solver parameters - INTEGER :: ITER, ITMAX,IERR,ITRACE, METHD,IPREC, ISTOPC,& - & IPARM(20), ML - REAL(KIND(1.D0)) :: ERR, EPS, RPARM(20) - - ! Other variables - INTEGER :: I,INFO - INTEGER :: INTERNAL, M,II - - ! Initialize BLACS - CALL BLACS_PINFO(IAM, NP) - CALL BLACS_GET(IZERO, IZERO, ICONTXT) - - ! Rectangular Grid, P x 1 - - CALL BLACS_GRIDINIT(ICONTXT, ORDER, NP, IONE) - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - ! - ! Get parameters - ! - CALL GET_PARMS(ICONTXT,CMETHD,PREC,AFMT,IDIM,ISTOPC,ITMAX,ITRACE,ML) - - ! - ! Allocate and fill in the coefficient matrix, RHS and initial guess - ! - - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - CALL CREATE_MATRIX(IDIM,A,B,X,DESC_A,PART_BLOCK,ICONTXT,AFMT) - T2 = MPI_WTIME() - T1 - - DIM=SIZE(A%ASPK) - - ALLOCATE(H%ASPK(DIM),H%IA1(DIM),H%IA2(DIM),H%PL(SIZE(A%PL)),& - & H%PL(SIZE(A%PL)),D(SIZE(A%PL)),& - & DESC_A_OUT%MATRIX_DATA(SIZE(DESC_A%MATRIX_DATA)),& - & DESC_A_OUT%HALO_INDEX(SIZE(DESC_A%HALO_INDEX)),& - & DESC_A_OUT%OVRLAP_INDEX(SIZE(DESC_A%OVRLAP_INDEX)),& - & DESC_A_OUT%OVRLAP_ELEM(SIZE(DESC_A%OVRLAP_ELEM)),& - & DESC_A_OUT%LOC_TO_GLOB(SIZE(DESC_A%LOC_TO_GLOB)),& - & DESC_A_OUT%GLOB_TO_LOC(SIZE(DESC_A%GLOB_TO_LOC)), WORK(1024)) - check_descr=15 -! work(5)=9 -!!$ WRITE(0,*)'CALLING VERIFY' -!!$ CALL F90_PSVERIFY(D,A,DESC_A,CHECK_DESCR,CONVERT_DESCR,H,& -!!$ & DESC_A_OUT,WORK) -!!$ WRITE(0,*)'VERIFY DONE',CONVERT_DESCR - - deallocate(work) - - CALL DGAMX2D(ICONTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1) - IF (IAM.EQ.0) Write(6,*) 'Matrix creation Time : ',T2 - - ! - ! Prepare the preconditioner. - ! - write(0,*)'PRECONDIZIONATORE=',prec - SELECT CASE (PREC) - CASE ('SCHW6') - IPREC = 6 - CASE ('SCHW5') - IPREC = 5 - CASE ('SCHW4') - IPREC = 4 - CASE ('SCHW3') - IPREC = 3 - CASE ('ILU') - IPREC = 2 - CASE ('DIAGSC') - IPREC = 1 - CASE ('NONE') - IPREC = 0 - CASE DEFAULT - WRITE(0,*) 'Unknown preconditioner' - CALL BLACS_ABORT(ICONTXT,-1) - END SELECT - pre%prec=iprec - pre%n_ovr=ml - pre%irenum=0 - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - CALL PRECONDITIONER(A,PRE,DESC_A,IRET) -!!$ CALL PRECONDITIONER(IPREC,A,L,U,D,DESC_A,IRET) - TPREC = MPI_WTIME()-T1 - - CALL DGAMX2D(icontxt,'A',' ',IONE, IONE,TPREC,IONE,t1,t1,-1,-1,-1) - - IF (IAM.EQ.0) WRITE(6,*) 'Preconditioner Time : ',TPREC - - IF (IRET.NE.0) THEN - WRITE(0,*) 'Error on preconditioner',IRET - CALL BLACS_ABORT(ICONTXT,-1) - STOP - END IF - - ! - ! Iterative method parameters - ! - call dcsprt90(80+myprow,a,head='% Local A') - - write(*,*) 'Calling Iterative method', size(b),ml - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - EPS = 1.D-9 - IF (CMETHD.EQ.'BICGSTAB') THEN - CALL F90_BICGSTAB(A,PRE,B,X,EPS,DESC_A,& - & ITMAX,ITER,ERR,IERR,ITRACE) -!!$ ELSE IF (CMETHD.EQ.'BICG') THEN -!!$ CALL F90_BICG(A,PRE,B,X,EPS,DESC_A,& -!!$ & ITMAX,ITER,ERR,IERR,ITRACE) - ELSE IF (CMETHD.EQ.'CGS') THEN - CALL F90_CGS(A,PRE,B,X,EPS,DESC_A,& - & ITMAX,ITER,ERR,IERR,ITRACE) - ELSE IF (CMETHD.EQ.'BICGSTABL') THEN - CALL F90_BICGSTABL(A,PRE,B,X,EPS,DESC_A,& - & ITMAX,ITER,ERR,IERR,ITRACE,ML) - ELSE - write(0,*) 'Unknown method ',cmethd - end IF - - CALL BLACS_BARRIER(ICONTXT,'All') - T2 = MPI_WTIME() - T1 - CALL DGAMX2D(ICONTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1) - - IF (IAM.EQ.0) THEN - WRITE(6,*) 'Time to Solve Matrix : ',T2 - WRITE(6,*) 'Time per iteration : ',T2/ITER - WRITE(6,*) 'Number of iterations : ',ITER - WRITE(6,*) 'Error on exit : ',ERR - WRITE(6,*) 'INFO on exit : ',IERR - END IF - - ! - ! Cleanup storage and exit - ! - CALL F90_PSDSFREE(B,DESC_A) - CALL F90_PSDSFREE(X,DESC_A) -!!$ CALL F90_PSDSFREE(D,DESC_A) - - CALL F90_PSSPFREE(A,DESC_A) -!!$ CALL F90_PSSPFREE(L,DESC_A) -!!$ CALL F90_PSSPFREE(U,DESC_A) - CALL F90_PSDSCFREE(DESC_A,info) - - CALL BLACS_GRIDEXIT(ICONTXT) - CALL BLACS_EXIT(0) - - STOP - -CONTAINS - ! - ! Get iteration parameters from the command line - ! - SUBROUTINE GET_PARMS(ICONTXT,CMETHD,PREC,AFMT,IDIM,ISTOPC,ITMAX,ITRACE,ML) - integer :: icontxt - Character :: CMETHD*10, PREC*10, AFMT*5 - Integer :: IDIM, IRET, ISTOPC,ITMAX,ITRACE,ML - Character*40 :: CHARBUF - INTEGER :: IARGC, NPROW, NPCOL, MYPROW, MYPCOL - EXTERNAL IARGC - INTEGER :: INTBUF(10), IP - - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - IF (MYPROW==0) THEN - READ(*,*) IP - IF (IP.GE.3) THEN - READ(*,*) CMETHD - READ(*,*) PREC - READ(*,*) AFMT - - ! Convert strings in array - DO I = 1, LEN(CMETHD) - INTBUF(I) = IACHAR(CMETHD(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - DO I = 1, LEN(PREC) - INTBUF(I) = IACHAR(PREC(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - DO I = 1, LEN(AFMT) - INTBUF(I) = IACHAR(AFMT(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - READ(*,*) IDIM - IF (IP.GE.4) THEN - READ(*,*) ISTOPC - ELSE - ISTOPC=1 - ENDIF - IF (IP.GE.5) THEN - READ(*,*) ITMAX - ELSE - ITMAX=500 - ENDIF - IF (IP.GE.6) THEN - READ(*,*) ITRACE - ELSE - ITRACE=-1 - ENDIF - IF (IP.GE.7) THEN - READ(*,*) ML - ELSE - ML=1 - ENDIF - ! Broadcast parameters to all processors - - INTBUF(1) = IDIM - INTBUF(2) = ISTOPC - INTBUF(3) = ITMAX - INTBUF(4) = ITRACE - INTBUF(5) = ML - CALL IGEBS2D(ICONTXT,'ALL',' ',5,1,INTBUF,5) - - WRITE(6,*)'Solving matrix: ELL1' - WRITE(6,*)'on grid',IDIM,'x',IDIM,'x',IDIM - WRITE(6,*)' with BLOCK data distribution, NP=',Np,& - & ' Preconditioner=',PREC,& - & ' Iterative methd=',CMETHD - ELSE - ! Wrong number of parameter, print an error message and exit - CALL PR_USAGE(0) - CALL BLACS_ABORT(ICONTXT,-1) - STOP 1 - ENDIF - ELSE - ! Receive Parameters - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 10 - CMETHD(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 10 - PREC(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 5 - AFMT(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',5,1,INTBUF,5,0,0) - IDIM = INTBUF(1) - ISTOPC = INTBUF(2) - ITMAX = INTBUF(3) - ITRACE = INTBUF(4) - ML = INTBUF(5) - END IF - RETURN - - END SUBROUTINE GET_PARMS - ! - ! Print an error message - ! - SUBROUTINE PR_USAGE(IOUT) - INTEGER :: IOUT - WRITE(IOUT,*)'Incorrect parameter(s) found' - WRITE(IOUT,*)' Usage: pde90 methd prec dim & - &[istop itmax itrace]' - WRITE(IOUT,*)' Where:' - WRITE(IOUT,*)' methd: CGSTAB TFQMR CGS' - WRITE(IOUT,*)' prec : ILU DIAGSC NONE' - WRITE(IOUT,*)' dim number of points along each axis' - WRITE(IOUT,*)' the size of the resulting linear ' - WRITE(IOUT,*)' system is dim**3' - WRITE(IOUT,*)' istop Stopping criterion 1, 2 or 3 [1] ' - WRITE(IOUT,*)' itmax Maximum number of iterations [500] ' - WRITE(IOUT,*)' itrace 0 (no tracing, default) or ' - WRITE(IOUT,*)' >= 0 do tracing every ITRACE' - WRITE(IOUT,*)' iterations ' - END SUBROUTINE PR_USAGE - -! -! Subroutine to allocate and fill in the coefficient matrix and -! the RHS. -! - SUBROUTINE CREATE_MATRIX(IDIM,A,B,T,DESC_A,PARTS,ICONTXT,AFMT) - ! - ! Discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) a1 d(u) a2 d(u) - ! - ------ - ------ - ----- - ------ + a4 u - ! dxdx dydy dx dy - ! - ! = 0 - ! - ! boundary condition: Dirichlet - ! 0< x,y<1 - ! - ! u(x,y)(2b1+2b2+a1+a2)+u(x-1,y)(-b1-a1)+u(x,y-1)(-b2-a2)+ - ! -u(x+1,y,z)b1-u(x,y+1,z)b2 - - USE TYPESP - USE TYPEDESC - USE F90TOOLS - Implicit None - INTEGER :: IDIM - integer, parameter :: nbmax=10 - Real(Kind(1.D0)),Pointer :: B(:),T(:) - Type (desc_type) :: DESC_A - Integer :: ICONTXT - INTERFACE - ! .....user passed subroutine..... - SUBROUTINE PARTS(GLOBAL_INDX,N,NP,PV,NV) - IMPLICIT NONE - INTEGER, INTENT(IN) :: GLOBAL_INDX, N, NP - INTEGER, INTENT(OUT) :: NV - INTEGER, INTENT(OUT) :: PV(*) - END SUBROUTINE PARTS - END INTERFACE ! Local variables - Type(D_SPMAT) :: A - Real(Kind(1.d0)) :: ZT(NBMAX),GLOB_X,GLOB_Y,GLOB_Z - Integer :: M,N,NNZ,GLOB_ROW,J - Type (D_SPMAT) :: ROW_MAT - Integer :: X,Y,Z,COUNTER,IA,I,INDX_OWNER - INTEGER :: NPROW,NPCOL,MYPROW,MYPCOL - Integer :: ELEMENT - INTEGER :: INFO, NV, INV - INTEGER, ALLOCATABLE :: PRV(:) - INTEGER, pointer :: ierrv(:) - Real(Kind(1.d0)), pointer :: DWORK(:) - INTEGER,POINTER :: IWORK(:) - character :: afmt*5 - ! deltah dimension of each grid cell - ! deltat discretization time - Real(Kind(1.D0)) :: DELTAH - Real(Kind(1.d0)),Parameter :: RHS=0.d0,ONE=1.d0,ZERO=0.d0 - Real(Kind(1.d0)) :: MPI_WTIME, T1, T2, T3, TINS - Real(Kind(1.d0)) :: a1, a2, a3, a4, b1, b2, b3 - external mpi_wtime,a1, a2, a3, a4, b1, b2, b3 - integer :: nb, ir1, ir2, ipr - logical :: own - ! common area - - - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - DELTAH = 1.D0/(IDIM-1) - - ! Initialize array descriptor and sparse matrix storage. Provide an - ! estimate of the number of non zeroes - CALL SETERR(2) - allocate(ierrv(6)) - - ierrv(:) = 0 - M = IDIM*IDIM - N = M - NNZ = ((N*6)/(NPROW*NPCOL)) - write(*,*) 'Size: n ',n - Call F90_PSDSCALL(N,N,PARTS,ICONTXT,IERRV,DESC_A) - write(*,*) 'Allocating A : nnz',nnz - Call F90_PSSPALL(A,IERRV,DESC_A,NNZ=NNZ) - ! Define RHS from boundary conditions; also build initial guess - write(*,*) 'Allocating B' - Call F90_PSDSALL(N,B,IERRV,DESC_A) - write(*,*) 'Allocating T' - Call F90_PSDSALL(N,T,IERRV,DESC_A) - - ! We build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. Might be extended to generate - ! a bunch of rows per call. - ! - ROW_MAT%DESCRA(1:1) = 'G' - ROW_MAT%FIDA = 'CSR' - write(*,*) 'Allocating ROW_MAT',20*nbmax - ALLOCATE(ROW_MAT%ASPK(20*nbmax),ROW_MAT%IA1(20*nbmax),& - &ROW_MAT%IA2(20*nbmax),PRV(NPROW),stat=info) - if (info.ne.0 ) then - write(*,*) 'Memory allocation error' - call blacs_abort(icontxt,-1) - endif - - TINS = 0.D0 - CALL BLACS_BARRIER(ICONTXT,'ALL') - T1 = MPI_WTIME() - - ! Loop over rows belonging to current process in a BLOCK - ! distribution. - - Z=0 - ROW_MAT%IA2(1)=1 - DO GLOB_ROW = 1, N - CALL PARTS(GLOB_ROW,N,NPROW,PRV,NV) - DO INV = 1, NV - INDX_OWNER = PRV(INV) - IF (INDX_OWNER == MYPROW) THEN - ! Local matrix pointer - ELEMENT=1 - ! Compute gridpoint Coordinates - IF (MOD(GLOB_ROW,(IDIM)).EQ.0) THEN - X = GLOB_ROW/(IDIM) - ELSE - X = GLOB_ROW/(IDIM)+1 - ENDIF - Y = GLOB_ROW-(X-1)*IDIM - ! GLOB_X, GLOB_Y, GLOB_X coordinates - GLOB_X=X*DELTAH - GLOB_Y=Y*DELTAH - GLOB_Z=Z*DELTAH - - - ! Check on boundary points -!!$ IF (X.EQ.1) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (Y.EQ.1) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (Z.EQ.1) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (X.EQ.IDIM) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (Y.EQ.IDIM) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (Z.EQ.IDIM) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE - zt(1) = 0.d0 - ! Internal point: build discretization - ! - ! Term depending on (x-1,y,z) - ! - if (x==1) then - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z)& - & -A1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*(-ROW_MAT%ASPK(ELEMENT)) - else - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z)& - & -A1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-2)*IDIM+(Y) - ELEMENT=ELEMENT+1 - endif - ! Term depending on (x,y-1,z) - if (y==1) then - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z)& - & -A2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)*(-ROW_MAT%ASPK(ELEMENT)) - else - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z)& - & -A2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y-1) - ELEMENT=ELEMENT+1 - endif - ! Term depending on (x,y,z-1) -!!$ if (z==1) then -!!$ ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z)& -!!$ & -A3(GLOB_X,GLOB_Y,GLOB_Z) -!!$ ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& -!!$ & DELTAH) -!!$ ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)*(-ROW_MAT%ASPK(ELEMENT)) -!!$ else -!!$ ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z)& -!!$ & -A3(GLOB_X,GLOB_Y,GLOB_Z) -!!$ ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& -!!$ & DELTAH) -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z-1) -!!$ ELEMENT=ELEMENT+1 -!!$ endif - ! Term depending on (x,y,z) - ROW_MAT%ASPK(ELEMENT)=2*B1(GLOB_X,GLOB_Y,GLOB_Z)& - & +2*B2(GLOB_X,GLOB_Y,GLOB_Z)& - & +A1(GLOB_X,GLOB_Y,GLOB_Z)& - & +A2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y) - ELEMENT=ELEMENT+1 - ! Term depending on (x,y,z+1) -!!$ if (z==idim) then -!!$ ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z) -!!$ ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& -!!$ & DELTAH) -!!$ ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)*(-ROW_MAT%ASPK(ELEMENT)) -!!$ else -!!$ ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z) -!!$ ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& -!!$ & DELTAH) -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z+1) -!!$ ELEMENT=ELEMENT+1 -!!$ endif - ! Term depending on (x,y+1,z) - if (y==idim) then - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)*(-ROW_MAT%ASPK(ELEMENT)) - else - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM+(Y+1) - ELEMENT=ELEMENT+1 - endif - ! Term depending on (x+1,y,z) - if (x= 0 do tracing every ITRACE' - WRITE(IOUT,*)' iterations ' - END SUBROUTINE PR_USAGE - -! -! Subroutine to allocate and fill in the coefficient matrix and -! the RHS. -! - SUBROUTINE CREATE_MATRIX(IDIM,A,B,T,DESC_A,PARTS,ICONTXT,AFMT) - ! - ! Discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) - ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u - ! dxdx dydy dzdz dx dy dz - ! - ! = 0 - ! - ! boundary condition: Dirichlet - ! 0< x,y,z<1 - ! - ! u(x,y,z)(2b1+2b2+2b3+a1+a2+a3)+u(x-1,y,z)(-b1-a1)+u(x,y-1,z)(-b2-a2)+ - ! + u(x,y,z-1)(-b3-a3)-u(x+1,y,z)b1-u(x,y+1,z)b2-u(x,y,z+1)b3 - - USE TYPESP - USE TYPEDESC - USE F90TOOLS - USE F90METHD - Implicit None - INTEGER :: IDIM - integer, parameter :: nbmax=10 - Real(Kind(1.D0)),Pointer :: B(:),T(:) - Type (desc_type) :: DESC_A - Integer :: ICONTXT - INTERFACE - ! .....user passed subroutine..... - SUBROUTINE PARTS(GLOBAL_INDX,N,NP,PV,NV) - IMPLICIT NONE - INTEGER, INTENT(IN) :: GLOBAL_INDX, N, NP - INTEGER, INTENT(OUT) :: NV - INTEGER, INTENT(OUT) :: PV(*) - END SUBROUTINE PARTS - END INTERFACE ! Local variables - Type(D_SPMAT) :: A - Real(Kind(1.d0)) :: ZT(NBMAX),GLOB_X,GLOB_Y,GLOB_Z - Integer :: M,N,NNZ,GLOB_ROW,J - Type (D_SPMAT) :: ROW_MAT - Integer :: X,Y,Z,COUNTER,IA,I,INDX_OWNER - INTEGER :: NPROW,NPCOL,MYPROW,MYPCOL - Integer :: ELEMENT - INTEGER :: INFO, NV, INV - INTEGER, ALLOCATABLE :: PRV(:) - INTEGER, pointer :: ierrv(:) - Real(Kind(1.d0)), pointer :: DWORK(:) - INTEGER,POINTER :: IWORK(:) - character :: afmt*5 - ! deltah dimension of each grid cell - ! deltat discretization time - Real(Kind(1.D0)) :: DELTAH - Real(Kind(1.d0)),Parameter :: RHS=0.d0,ONE=1.d0,ZERO=0.d0 - Real(Kind(1.d0)) :: MPI_WTIME, T1, T2, T3, TINS - Real(Kind(1.d0)) :: a1, a2, a3, a4, b1, b2, b3 - external mpi_wtime,a1, a2, a3, a4, b1, b2, b3 - integer :: nb, ir1, ir2, ipr - logical :: own - ! common area - - - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - DELTAH = 1.D0/(IDIM-1) - - ! Initialize array descriptor and sparse matrix storage. Provide an - ! estimate of the number of non zeroes - CALL SETERR(2) - allocate(ierrv(6)) - - ierrv(:) = 0 - M = IDIM*IDIM*IDIM - N = M - NNZ = ((N*9)/(NPROW*NPCOL)) - write(*,*) 'Size: n ',n - Call F90_PSDSCALL(N,N,PARTS,ICONTXT,IERRV,DESC_A) - write(*,*) 'Allocating A : nnz',nnz - Call F90_PSSPALL(A,IERRV,DESC_A,NNZ=NNZ) - ! Define RHS from boundary conditions; also build initial guess - write(*,*) 'Allocating B' - Call F90_PSDSALL(N,B,IERRV,DESC_A) - write(*,*) 'Allocating T' - Call F90_PSDSALL(N,T,IERRV,DESC_A) - - ! We build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. Might be extended to generate - ! a bunch of rows per call. - ! - ROW_MAT%DESCRA(1:1) = 'G' - ROW_MAT%FIDA = 'CSR' - write(*,*) 'Allocating ROW_MAT',20*nbmax - ALLOCATE(ROW_MAT%ASPK(20*nbmax),ROW_MAT%IA1(20*nbmax),& - &ROW_MAT%IA2(20*nbmax),PRV(NPROW),stat=info) - if (info.ne.0 ) then - write(*,*) 'Memory allocation error' - call blacs_abort(icontxt,-1) - endif - - TINS = 0.D0 - CALL BLACS_BARRIER(ICONTXT,'ALL') - T1 = MPI_WTIME() - - ! Loop over rows belonging to current process in a BLOCK - ! distribution. - - ROW_MAT%IA2(1)=1 - DO GLOB_ROW = 1, N - CALL PARTS(GLOB_ROW,N,NPROW,PRV,NV) - DO INV = 1, NV - INDX_OWNER = PRV(INV) - IF (INDX_OWNER == MYPROW) THEN - ! Local matrix pointer - ELEMENT=1 - ! Compute gridpoint Coordinates - IF (MOD(GLOB_ROW,(IDIM*IDIM)).EQ.0) THEN - X = GLOB_ROW/(IDIM*IDIM) - ELSE - X = GLOB_ROW/(IDIM*IDIM)+1 - ENDIF - IF (MOD((GLOB_ROW-(X-1)*IDIM*IDIM),IDIM).EQ.0) THEN - Y = (GLOB_ROW-(X-1)*IDIM*IDIM)/IDIM - ELSE - Y = (GLOB_ROW-(X-1)*IDIM*IDIM)/IDIM+1 - ENDIF - Z = GLOB_ROW-(X-1)*IDIM*IDIM-(Y-1)*IDIM - ! GLOB_X, GLOB_Y, GLOB_X coordinates - GLOB_X=X*DELTAH - GLOB_Y=Y*DELTAH - GLOB_Z=Z*DELTAH - - ! Check on boundary points - IF (X.EQ.1) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ELSE IF (Y.EQ.1) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ELSE IF (Z.EQ.1) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ELSE IF (X.EQ.IDIM) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ELSE IF (Y.EQ.IDIM) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ELSE IF (Z.EQ.IDIM) THEN - ROW_MAT%ASPK(ELEMENT)=ONE - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ELSE - ! Internal point: build discretization - ! - ! Term depending on (x-1,y,z) - ! - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z)& - & -A1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-2)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ! Term depending on (x,y-1,z) - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z)& - & -A2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-2)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ! Term depending on (x,y,z-1) - ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z)& - & -A3(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z-1) - ELEMENT=ELEMENT+1 - ! Term depending on (x,y,z) - ROW_MAT%ASPK(ELEMENT)=2*B1(GLOB_X,GLOB_Y,GLOB_Z)& - & +2*B2(GLOB_X,GLOB_Y,GLOB_Z)& - & +2*B3(GLOB_X,GLOB_Y,GLOB_Z)& - & +A1(GLOB_X,GLOB_Y,GLOB_Z)& - & +A2(GLOB_X,GLOB_Y,GLOB_Z)& - & +A3(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ! Term depending on (x,y,z+1) - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z+1) - ELEMENT=ELEMENT+1 - ! Term depending on (x,y+1,z) - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ! Term depending on (x+1,y,z) - ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ENDIF - ROW_MAT%M=1 - ROW_MAT%K=N - ROW_MAT%IA2(2)=ELEMENT - ! IA== GLOBAL ROW INDEX - IA=GLOB_ROW -!!$ IA=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ write(0,*) 'Inserting row ',ia,' On proc',myprow - T3 = MPI_WTIME() - CALL F90_PSSPINS(A,IA,1,ROW_MAT,IERRV,DESC_A) - if (ierrv(1).ne.0) then - write(0,*) 'On row ',ia,' IERRV:',ierrv(:) - endif - TINS = TINS + (MPI_WTIME()-T3) - ! Build RHS - IF (X==1) THEN - GLOB_Y=(Y-IDIM/2)*DELTAH - GLOB_Z=(Z-IDIM/2)*DELTAH - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2) - ELSE IF ((Y==1).OR.(Y==IDIM).OR.(Z==1).OR.(Z==IDIM)) THEN - GLOB_X=3*(X-1)*DELTAH - GLOB_Y=(Y-IDIM/2)*DELTAH - GLOB_Z=(Z-IDIM/2)*DELTAH - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X) - ELSE - ZT(1) = 0.D0 - ENDIF - CALL F90_PSDSINS(1,B,IA,ZT(1:1),IERRV,DESC_A) - ZT(1)=0.D0 - CALL F90_PSDSINS(1,T,IA,ZT(1:1),IERRV,DESC_A) - END IF - END DO - END DO - - CALL BLACS_BARRIER(ICONTXT,'ALL') - T2 = MPI_WTIME() - - WRITE(*,*) ' pspins time',TINS - WRITE(*,*) ' Insert time',(T2-T1) - - DEALLOCATE(ROW_MAT%ASPK,ROW_MAT%IA1,ROW_MAT%IA2) - - write(*,*) 'Calling SPASB' - CALL BLACS_BARRIER(ICONTXT,'ALL') - T1 = MPI_WTIME() - - CALL F90_PSSPASB(A,IERRV,DESC_A,AFMT=AFMT,DUP=2) - - CALL BLACS_BARRIER(ICONTXT,'ALL') - T2 = MPI_WTIME() - - WRITE(0,*) ' Assembly time',(T2-T1),' ',a%fida(1:4) - - CALL F90_PSDSASB(B,IERRV,DESC_A) - CALL F90_PSDSASB(T,IERRV,DESC_A) - IF (MYPROW.EQ.0) THEN - WRITE(0,*) ' End CREATE_MATRIX' - ENDIF - RETURN - - END SUBROUTINE CREATE_MATRIX -END PROGRAM PDE90 -! -! Functions parametrizing the differential equation -! -FUNCTION A1(X,Y,Z) - REAL(KIND(1.D0)) :: A1 - REAL(KIND(1.D0)) :: X,Y,Z - A1=1.D0 -END FUNCTION A1 -FUNCTION A2(X,Y,Z) - REAL(KIND(1.D0)) :: A2 - REAL(KIND(1.D0)) :: X,Y,Z - A2=2.D1*Y -END FUNCTION A2 -FUNCTION A3(X,Y,Z) - REAL(KIND(1.D0)) :: A3 - REAL(KIND(1.D0)) :: X,Y,Z - A3=1.D0 -END FUNCTION A3 -FUNCTION A4(X,Y,Z) - REAL(KIND(1.D0)) :: A4 - REAL(KIND(1.D0)) :: X,Y,Z - A4=1.D0 -END FUNCTION A4 -FUNCTION B1(X,Y,Z) - REAL(KIND(1.D0)) :: B1 - REAL(KIND(1.D0)) :: X,Y,Z - B1=1.D0 -END FUNCTION B1 -FUNCTION B2(X,Y,Z) - REAL(KIND(1.D0)) :: B2 - REAL(KIND(1.D0)) :: X,Y,Z - B2=1.D0 -END FUNCTION B2 -FUNCTION B3(X,Y,Z) - REAL(KIND(1.D0)) :: B3 - REAL(KIND(1.D0)) :: X,Y,Z - B3=1.D0 -END FUNCTION B3 - - diff --git a/test/pargen/ppde90s.f90 b/test/pargen/ppde90s.f90 deleted file mode 100644 index eb1e3d22..00000000 --- a/test/pargen/ppde90s.f90 +++ /dev/null @@ -1,707 +0,0 @@ -! -! This sample program shows how to build and solve a sparse linear -! -! The program solves a linear system based on the partial differential -! equation -! -! -! -! the equation generated is: -! b1 d d (u) b2 d d (u) a1 d (u)) a2 d (u))) -! - ------ - ------ + ----- + ------ + a3 u = 0 -! dx dx dy dy dx dy -! -! -! with Dirichlet boundary conditions on the unit cube -! -! 0<=x,y,z<=1 -! -! The equation is discretized with finite differences and uniform stepsize; -! the resulting discrete equation is -! -! ( u(x,y,z)(2b1+2b2+a1+a2)+u(x-1,y)(-b1-a1)+u(x,y-1)(-b2-a2)+ -! -u(x+1,y)b1-u(x,y+1)b2)*(1/h**2) -! -! Example taken from: C.T.Kelley -! Iterative Methods for Linear and Nonlinear Equations -! SIAM 1995 -! -! -! In this sample program the index space of the discretized -! computational domain is first numbered sequentially in a standard way, -! then the corresponding vector is distributed according to an HPF BLOCK -! distribution directive. -! -! Boundary conditions are set in a very simple way, by adding -! equations of the form -! -! u(x,y) = rhs(x,y) -! -Program PDE90 - USE F90SPARSE - Implicit none - - interface - !.....user passed subroutine..... - subroutine part_block(glob_index,n,np,pv,nv) - INTEGER, INTENT(IN) :: GLOB_INDEX, N, NP - INTEGER, INTENT(OUT) :: NV - INTEGER, INTENT(OUT) :: PV(*) - end subroutine part_block - end interface - ! input parameters - Character :: CMETHD*10, PREC*10, AFMT*5 - Integer :: IDIM, IRET - - ! Miscellaneous - Integer, Parameter :: IZERO=0, IONE=1 - Character, PARAMETER :: ORDER='R' - INTEGER :: IARGC,CONVERT_DESCR,dim, CHECK_DESCR - REAL(KIND(1.D0)), PARAMETER :: DZERO = 0.D0, ONE = 1.D0 - REAL(KIND(1.D0)) :: MPI_WTIME, T1, T2, TPREC, TSOLVE, T3, T4 - EXTERNAL MPI_WTIME - - ! Sparse Matrix and preconditioner - TYPE(D_SPMAT) :: A, L, U, H - TYPE(D_PREC) :: PRE - ! Descriptor - TYPE(desc_type) :: DESC_A, DESC_A_OUT - ! Dense Matrices - REAL(KIND(1.d0)), POINTER :: B(:), X(:), D(:),LD(:) - INTEGER, pointer :: WORK(:) - ! BLACS parameters - INTEGER :: nprow, npcol, icontxt, iam, np, myprow, mypcol - - ! Solver parameters - INTEGER :: ITER, ITMAX,IERR,ITRACE, METHD,IPREC, ISTOPC,& - & IPARM(20), ML - REAL(KIND(1.D0)) :: ERR, EPS, RPARM(20) - - ! Other variables - INTEGER :: I,INFO - INTEGER :: INTERNAL, M,II - - ! Initialize BLACS - CALL BLACS_PINFO(IAM, NP) - CALL BLACS_GET(IZERO, IZERO, ICONTXT) - - ! Rectangular Grid, P x 1 - - CALL BLACS_GRIDINIT(ICONTXT, ORDER, NP, IONE) - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - ! - ! Get parameters - ! - CALL GET_PARMS(ICONTXT,CMETHD,PREC,AFMT,IDIM,ISTOPC,ITMAX,ITRACE,ML) - - ! - ! Allocate and fill in the coefficient matrix, RHS and initial guess - ! - - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - CALL CREATE_MATRIX(IDIM,A,B,X,DESC_A,PART_BLOCK,ICONTXT,AFMT) - T2 = MPI_WTIME() - T1 - - DIM=SIZE(A%ASPK) - - ALLOCATE(H%ASPK(DIM),H%IA1(DIM),H%IA2(DIM),H%PL(SIZE(A%PL)),& - & H%PL(SIZE(A%PL)),D(SIZE(A%PL)),& - & DESC_A_OUT%MATRIX_DATA(SIZE(DESC_A%MATRIX_DATA)),& - & DESC_A_OUT%HALO_INDEX(SIZE(DESC_A%HALO_INDEX)),& - & DESC_A_OUT%OVRLAP_INDEX(SIZE(DESC_A%OVRLAP_INDEX)),& - & DESC_A_OUT%OVRLAP_ELEM(SIZE(DESC_A%OVRLAP_ELEM)),& - & DESC_A_OUT%LOC_TO_GLOB(SIZE(DESC_A%LOC_TO_GLOB)),& - & DESC_A_OUT%GLOB_TO_LOC(SIZE(DESC_A%GLOB_TO_LOC)), WORK(1024)) - check_descr=15 -! work(5)=9 -!!$ WRITE(0,*)'CALLING VERIFY' -!!$ CALL F90_PSVERIFY(D,A,DESC_A,CHECK_DESCR,CONVERT_DESCR,H,& -!!$ & DESC_A_OUT,WORK) -!!$ WRITE(0,*)'VERIFY DONE',CONVERT_DESCR - - deallocate(work) - - CALL DGAMX2D(ICONTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1) - IF (IAM.EQ.0) Write(6,*) 'Matrix creation Time : ',T2 - - ! - ! Prepare the preconditioner. - ! - write(0,*)'PRECONDIZIONATORE=',prec - SELECT CASE (PREC) - CASE ('SCHW6') - IPREC = 6 - CASE ('SCHW5') - IPREC = 5 - CASE ('SCHW4') - IPREC = 4 - CASE ('SCHW3') - IPREC = 3 - CASE ('ILU') - IPREC = 2 - CASE ('DIAGSC') - IPREC = 1 - CASE ('NONE') - IPREC = 0 - CASE DEFAULT - WRITE(0,*) 'Unknown preconditioner' - CALL BLACS_ABORT(ICONTXT,-1) - END SELECT - pre%prec=iprec - pre%n_ovr=ml - pre%irenum=0 - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - CALL PRECONDITIONER(A,PRE,DESC_A,IRET) -!!$ CALL PRECONDITIONER(IPREC,A,L,U,D,DESC_A,IRET) - TPREC = MPI_WTIME()-T1 - - CALL DGAMX2D(icontxt,'A',' ',IONE, IONE,TPREC,IONE,t1,t1,-1,-1,-1) - - IF (IAM.EQ.0) WRITE(6,*) 'Preconditioner Time : ',TPREC - - IF (IRET.NE.0) THEN - WRITE(0,*) 'Error on preconditioner',IRET - CALL BLACS_ABORT(ICONTXT,-1) - STOP - END IF - - ! - ! Iterative method parameters - ! - call dcsprt90(80+myprow,a,head='% Local A') - - write(*,*) 'Calling Iterative method', size(b),ml - CALL BLACS_BARRIER(ICONTXT,'All') - T1 = MPI_WTIME() - EPS = 1.D-9 - IF (CMETHD.EQ.'BICGSTAB') THEN - CALL F90_BICGSTAB(A,PRE,B,X,EPS,DESC_A,& - & ITMAX,ITER,ERR,IERR,ITRACE) -!!$ ELSE IF (CMETHD.EQ.'BICG') THEN -!!$ CALL F90_BICG(A,PRE,B,X,EPS,DESC_A,& -!!$ & ITMAX,ITER,ERR,IERR,ITRACE) - ELSE IF (CMETHD.EQ.'CGS') THEN - CALL F90_CGS(A,PRE,B,X,EPS,DESC_A,& - & ITMAX,ITER,ERR,IERR,ITRACE) - ELSE IF (CMETHD.EQ.'BICGSTABL') THEN - CALL F90_BICGSTABL(A,PRE,B,X,EPS,DESC_A,& - & ITMAX,ITER,ERR,IERR,ITRACE,ML) - ELSE - write(0,*) 'Unknown method ',cmethd - end IF - - CALL BLACS_BARRIER(ICONTXT,'All') - T2 = MPI_WTIME() - T1 - CALL DGAMX2D(ICONTXT,'A',' ',IONE, IONE,T2,IONE,T1,T1,-1,-1,-1) - - IF (IAM.EQ.0) THEN - WRITE(6,*) 'Time to Solve Matrix : ',T2 - WRITE(6,*) 'Time per iteration : ',T2/ITER - WRITE(6,*) 'Number of iterations : ',ITER - WRITE(6,*) 'Error on exit : ',ERR - WRITE(6,*) 'INFO on exit : ',IERR - END IF - - ! - ! Cleanup storage and exit - ! - CALL F90_PSDSFREE(B,DESC_A) - CALL F90_PSDSFREE(X,DESC_A) -!!$ CALL F90_PSDSFREE(D,DESC_A) - - CALL F90_PSSPFREE(A,DESC_A) -!!$ CALL F90_PSSPFREE(L,DESC_A) -!!$ CALL F90_PSSPFREE(U,DESC_A) - CALL F90_PSDSCFREE(DESC_A,info) - - CALL BLACS_GRIDEXIT(ICONTXT) - CALL BLACS_EXIT(0) - - STOP - -CONTAINS - ! - ! Get iteration parameters from the command line - ! - SUBROUTINE GET_PARMS(ICONTXT,CMETHD,PREC,AFMT,IDIM,ISTOPC,ITMAX,ITRACE,ML) - integer :: icontxt - Character :: CMETHD*10, PREC*10, AFMT*5 - Integer :: IDIM, IRET, ISTOPC,ITMAX,ITRACE,ML - Character*40 :: CHARBUF - INTEGER :: IARGC, NPROW, NPCOL, MYPROW, MYPCOL - EXTERNAL IARGC - INTEGER :: INTBUF(10), IP - - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - IF (MYPROW==0) THEN - READ(*,*) IP - IF (IP.GE.3) THEN - READ(*,*) CMETHD - READ(*,*) PREC - READ(*,*) AFMT - - ! Convert strings in array - DO I = 1, LEN(CMETHD) - INTBUF(I) = IACHAR(CMETHD(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - DO I = 1, LEN(PREC) - INTBUF(I) = IACHAR(PREC(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - DO I = 1, LEN(AFMT) - INTBUF(I) = IACHAR(AFMT(I:I)) - END DO - ! Broadcast parameters to all processors - CALL IGEBS2D(ICONTXT,'ALL',' ',10,1,INTBUF,10) - - READ(*,*) IDIM - IF (IP.GE.4) THEN - READ(*,*) ISTOPC - ELSE - ISTOPC=1 - ENDIF - IF (IP.GE.5) THEN - READ(*,*) ITMAX - ELSE - ITMAX=500 - ENDIF - IF (IP.GE.6) THEN - READ(*,*) ITRACE - ELSE - ITRACE=-1 - ENDIF - IF (IP.GE.7) THEN - READ(*,*) ML - ELSE - ML=1 - ENDIF - ! Broadcast parameters to all processors - - INTBUF(1) = IDIM - INTBUF(2) = ISTOPC - INTBUF(3) = ITMAX - INTBUF(4) = ITRACE - INTBUF(5) = ML - CALL IGEBS2D(ICONTXT,'ALL',' ',5,1,INTBUF,5) - - WRITE(6,*)'Solving matrix: ELL1' - WRITE(6,*)'on grid',IDIM,'x',IDIM,'x',IDIM - WRITE(6,*)' with BLOCK data distribution, NP=',Np,& - & ' Preconditioner=',PREC,& - & ' Iterative methd=',CMETHD - ELSE - ! Wrong number of parameter, print an error message and exit - CALL PR_USAGE(0) - CALL BLACS_ABORT(ICONTXT,-1) - STOP 1 - ENDIF - ELSE - ! Receive Parameters - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 10 - CMETHD(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 10 - PREC(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',10,1,INTBUF,10,0,0) - DO I = 1, 5 - AFMT(I:I) = ACHAR(INTBUF(I)) - END DO - CALL IGEBR2D(ICONTXT,'ALL',' ',5,1,INTBUF,5,0,0) - IDIM = INTBUF(1) - ISTOPC = INTBUF(2) - ITMAX = INTBUF(3) - ITRACE = INTBUF(4) - ML = INTBUF(5) - END IF - RETURN - - END SUBROUTINE GET_PARMS - ! - ! Print an error message - ! - SUBROUTINE PR_USAGE(IOUT) - INTEGER :: IOUT - WRITE(IOUT,*)'Incorrect parameter(s) found' - WRITE(IOUT,*)' Usage: pde90 methd prec dim & - &[istop itmax itrace]' - WRITE(IOUT,*)' Where:' - WRITE(IOUT,*)' methd: CGSTAB TFQMR CGS' - WRITE(IOUT,*)' prec : ILU DIAGSC NONE' - WRITE(IOUT,*)' dim number of points along each axis' - WRITE(IOUT,*)' the size of the resulting linear ' - WRITE(IOUT,*)' system is dim**3' - WRITE(IOUT,*)' istop Stopping criterion 1, 2 or 3 [1] ' - WRITE(IOUT,*)' itmax Maximum number of iterations [500] ' - WRITE(IOUT,*)' itrace 0 (no tracing, default) or ' - WRITE(IOUT,*)' >= 0 do tracing every ITRACE' - WRITE(IOUT,*)' iterations ' - END SUBROUTINE PR_USAGE - -! -! Subroutine to allocate and fill in the coefficient matrix and -! the RHS. -! - SUBROUTINE CREATE_MATRIX(IDIM,A,B,T,DESC_A,PARTS,ICONTXT,AFMT) - ! - ! Discretize the partial diferential equation - ! - ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) - ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u - ! dxdx dydy dzdz dx dy dz - ! - ! = 0 - ! - ! boundary condition: Dirichlet - ! 0< x,y,z<1 - ! - ! u(x,y,z)(2b1+2b2+2b3+a1+a2+a3)+u(x-1,y,z)(-b1-a1)+u(x,y-1,z)(-b2-a2)+ - ! + u(x,y,z-1)(-b3-a3)-u(x+1,y,z)b1-u(x,y+1,z)b2-u(x,y,z+1)b3 - - USE TYPESP - USE TYPEDESC - USE F90TOOLS - USE F90METHD - Implicit None - INTEGER :: IDIM - integer, parameter :: nbmax=10 - Real(Kind(1.D0)),Pointer :: B(:),T(:) - Type (desc_type) :: DESC_A - Integer :: ICONTXT - INTERFACE - ! .....user passed subroutine..... - SUBROUTINE PARTS(GLOBAL_INDX,N,NP,PV,NV) - IMPLICIT NONE - INTEGER, INTENT(IN) :: GLOBAL_INDX, N, NP - INTEGER, INTENT(OUT) :: NV - INTEGER, INTENT(OUT) :: PV(*) - END SUBROUTINE PARTS - END INTERFACE ! Local variables - Type(D_SPMAT) :: A - Real(Kind(1.d0)) :: ZT(NBMAX),GLOB_X,GLOB_Y,GLOB_Z - Integer :: M,N,NNZ,GLOB_ROW,J - Type (D_SPMAT) :: ROW_MAT - Integer :: X,Y,Z,COUNTER,IA,I,INDX_OWNER - INTEGER :: NPROW,NPCOL,MYPROW,MYPCOL - Integer :: ELEMENT - INTEGER :: INFO, NV, INV - INTEGER, ALLOCATABLE :: PRV(:) - INTEGER, pointer :: ierrv(:) - Real(Kind(1.d0)), pointer :: DWORK(:) - INTEGER,POINTER :: IWORK(:) - character :: afmt*5 - ! deltah dimension of each grid cell - ! deltat discretization time - Real(Kind(1.D0)) :: DELTAH - Real(Kind(1.d0)),Parameter :: RHS=0.d0,ONE=1.d0,ZERO=0.d0 - Real(Kind(1.d0)) :: MPI_WTIME, T1, T2, T3, TINS - Real(Kind(1.d0)) :: a1, a2, a3, a4, b1, b2, b3 - external mpi_wtime,a1, a2, a3, a4, b1, b2, b3 - integer :: nb, ir1, ir2, ipr - logical :: own - ! common area - - - CALL BLACS_GRIDINFO(ICONTXT, NPROW, NPCOL, MYPROW, MYPCOL) - - DELTAH = 1.D0/(IDIM-1) - - ! Initialize array descriptor and sparse matrix storage. Provide an - ! estimate of the number of non zeroes - CALL SETERR(2) - allocate(ierrv(6)) - - ierrv(:) = 0 - M = IDIM*IDIM*IDIM - N = M - NNZ = ((N*9)/(NPROW*NPCOL)) - write(*,*) 'Size: n ',n - Call F90_PSDSCALL(N,N,PARTS,ICONTXT,IERRV,DESC_A) - write(*,*) 'Allocating A : nnz',nnz - Call F90_PSSPALL(A,IERRV,DESC_A,NNZ=NNZ) - ! Define RHS from boundary conditions; also build initial guess - write(*,*) 'Allocating B' - Call F90_PSDSALL(N,B,IERRV,DESC_A) - write(*,*) 'Allocating T' - Call F90_PSDSALL(N,T,IERRV,DESC_A) - - ! We build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. Might be extended to generate - ! a bunch of rows per call. - ! - ROW_MAT%DESCRA(1:1) = 'G' - ROW_MAT%FIDA = 'CSR' - write(*,*) 'Allocating ROW_MAT',20*nbmax - ALLOCATE(ROW_MAT%ASPK(20*nbmax),ROW_MAT%IA1(20*nbmax),& - &ROW_MAT%IA2(20*nbmax),PRV(NPROW),stat=info) - if (info.ne.0 ) then - write(*,*) 'Memory allocation error' - call blacs_abort(icontxt,-1) - endif - - TINS = 0.D0 - CALL BLACS_BARRIER(ICONTXT,'ALL') - T1 = MPI_WTIME() - - ! Loop over rows belonging to current process in a BLOCK - ! distribution. - - ROW_MAT%IA2(1)=1 - DO GLOB_ROW = 1, N - CALL PARTS(GLOB_ROW,N,NPROW,PRV,NV) - DO INV = 1, NV - INDX_OWNER = PRV(INV) - IF (INDX_OWNER == MYPROW) THEN - ! Local matrix pointer - ELEMENT=1 - ! Compute gridpoint Coordinates - IF (MOD(GLOB_ROW,(IDIM*IDIM)).EQ.0) THEN - X = GLOB_ROW/(IDIM*IDIM) - ELSE - X = GLOB_ROW/(IDIM*IDIM)+1 - ENDIF - IF (MOD((GLOB_ROW-(X-1)*IDIM*IDIM),IDIM).EQ.0) THEN - Y = (GLOB_ROW-(X-1)*IDIM*IDIM)/IDIM - ELSE - Y = (GLOB_ROW-(X-1)*IDIM*IDIM)/IDIM+1 - ENDIF - Z = GLOB_ROW-(X-1)*IDIM*IDIM-(Y-1)*IDIM - ! GLOB_X, GLOB_Y, GLOB_X coordinates - GLOB_X=X*DELTAH - GLOB_Y=Y*DELTAH - GLOB_Z=Z*DELTAH - - - ! Check on boundary points -!!$ IF (X.EQ.1) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (Y.EQ.1) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (Z.EQ.1) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (X.EQ.IDIM) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (Y.EQ.IDIM) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE IF (Z.EQ.IDIM) THEN -!!$ ROW_MAT%ASPK(ELEMENT)=ONE -!!$ ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) -!!$ ELEMENT=ELEMENT+1 -!!$ ELSE - zt(1) = 0.d0 - ! Internal point: build discretization - ! - ! Term depending on (x-1,y,z) - ! - if (x==1) then - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z)& - & -A1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*(-ROW_MAT%ASPK(ELEMENT)) - else - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z)& - & -A1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-2)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - endif - ! Term depending on (x,y-1,z) - if (y==1) then - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z)& - & -A2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)*(-ROW_MAT%ASPK(ELEMENT)) - else - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z)& - & -A2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-2)*IDIM+(Z) - ELEMENT=ELEMENT+1 - endif - ! Term depending on (x,y,z-1) - if (z==1) then - ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z)& - & -A3(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)*(-ROW_MAT%ASPK(ELEMENT)) - else - ROW_MAT%ASPK(ELEMENT)=-B3(GLOB_X,GLOB_Y,GLOB_Z)& - & -A3(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z-1) - ELEMENT=ELEMENT+1 - endif - ! Term depending on (x,y,z) - ROW_MAT%ASPK(ELEMENT)=2*B1(GLOB_X,GLOB_Y,GLOB_Z)& - & +2*B2(GLOB_X,GLOB_Y,GLOB_Z)& - & +2*B3(GLOB_X,GLOB_Y,GLOB_Z)& - & +A1(GLOB_X,GLOB_Y,GLOB_Z)& - & +A2(GLOB_X,GLOB_Y,GLOB_Z)& - & +A3(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z) - ELEMENT=ELEMENT+1 - ! Term depending on (x,y,z+1) - if (z==idim) then - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)*(-ROW_MAT%ASPK(ELEMENT)) - else - ROW_MAT%ASPK(ELEMENT)=-B1(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y-1)*IDIM+(Z+1) - ELEMENT=ELEMENT+1 - endif - ! Term depending on (x,y+1,z) - if (y==idim) then - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ZT(1) = EXP(-GLOB_Y**2-GLOB_Z**2)*EXP(-GLOB_X)*(-ROW_MAT%ASPK(ELEMENT)) - else - ROW_MAT%ASPK(ELEMENT)=-B2(GLOB_X,GLOB_Y,GLOB_Z) - ROW_MAT%ASPK(ELEMENT) = ROW_MAT%ASPK(ELEMENT)/(DELTAH*& - & DELTAH) - ROW_MAT%IA1(ELEMENT)=(X-1)*IDIM*IDIM+(Y)*IDIM+(Z) - ELEMENT=ELEMENT+1 - endif - ! Term depending on (x+1,y,z) - if (x