cleanup
Alfredo Buttari
19 years ago
parent
4745a96f43
commit
8919dacf49
@ -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
|
|
||||||
|
|
||||||
|
|
||||||
@ -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<idim) then
|
|
||||||
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
|
|
||||||
!!$ 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
|
|
||||||
|
|
||||||
|
|
||||||
@ -1,682 +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 TYPESP
|
|
||||||
USE TYPEDESC
|
|
||||||
USE TYPEPREC
|
|
||||||
USE F90TOOLS
|
|
||||||
USE F90METHD
|
|
||||||
USE F90PREC
|
|
||||||
use mpi
|
|
||||||
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)) :: T1, T2, TPREC, TSOLVE, T3, T4
|
|
||||||
!!$ EXTERNAL MPI_WTIME
|
|
||||||
integer mpe_log_get_event_number,mpe_Describe_state,mpe_log_event
|
|
||||||
|
|
||||||
! Sparse Matrix and preconditioner
|
|
||||||
TYPE(D_SPMAT) :: A, L, U, H
|
|
||||||
TYPE(PREC_DATA) :: 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,iprecb,iprece,islvb,islve
|
|
||||||
INTEGER :: INTERNAL, M,II
|
|
||||||
|
|
||||||
! Initialize BLACS
|
|
||||||
CALL BLACS_PINFO(IAM, NP)
|
|
||||||
CALL BLACS_GET(IZERO, IZERO, ICONTXT)
|
|
||||||
iprecb = mpe_log_get_event_number()
|
|
||||||
iprece = mpe_log_get_event_number()
|
|
||||||
islvb = mpe_log_get_event_number()
|
|
||||||
islve = mpe_log_get_event_number()
|
|
||||||
if (iam==0) then
|
|
||||||
info = mpe_describe_state(iprecb,iprece,"Preconditioner","OrangeRed")
|
|
||||||
info = mpe_describe_state(islvb,islve,"Solver","DarkGreen")
|
|
||||||
endif
|
|
||||||
|
|
||||||
|
|
||||||
! 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
|
|
||||||
CALL BLACS_BARRIER(ICONTXT,'All')
|
|
||||||
T1 = MPI_WTIME()
|
|
||||||
info = MPE_Log_event( iprecb, 0, "start Precond" )
|
|
||||||
|
|
||||||
CALL PRECONDITIONER(A,PRE,DESC_A,IRET)
|
|
||||||
!!$ CALL PRECONDITIONER(IPREC,A,L,U,D,DESC_A,IRET)
|
|
||||||
info = MPE_Log_event( iprece, 0, "end Precond" )
|
|
||||||
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
|
|
||||||
info = MPE_Log_event( islvb, 0, "start Solver" )
|
|
||||||
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
|
|
||||||
info = MPE_Log_event( islve, 0, "end Solver" )
|
|
||||||
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
|
|
||||||
! 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
|
|
||||||
|
|
||||||
|
|
||||||
@ -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<idim) then
|
|
||||||
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
|
|
||||||
!!$ 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
|
|
||||||
|
|
||||||
|
|
||||||
Loading…
Reference in New Issue