PROGRAM petscitest IMPLICIT NONE !To run code: mpiexec -n 1 ./petsctest -mat_view_info ! The following include statements are required for KSP Fortran programs: ! petscsys.h - base PETSc routines ! petscvec.h - vectors ! petscmat.h - matrices ! petscpc.h - preconditioners ! petscksp.h - Krylov subspace methods #include #include #include #include #include KSP :: cgsolve KSPConvergedReason :: reason Mat :: A Vec :: x,y,b PetscErrorCode :: ierr PetscMPIInt :: rank,size PetscInt :: N2 PetscScalar :: value PetscInt :: nr,np,nt,npts,num_steps,ndiags,ai,i,j,k,M,N PetscInt :: istart,iend integer, dimension(:), allocatable :: veci real*8, dimension(:), allocatable :: xvec real*8, dimension(:), allocatable :: yvec real*8, dimension(:), allocatable :: CSR_A integer, dimension(:), allocatable :: CSR_AI integer, dimension(:), allocatable :: CSR_AJ real ttmp1,ttmp2,flops,CSRtime,PETSctime real elapsed(2) real etime call PetscInitialize(PETSC_NULL_CHARACTER,ierr) call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr) npts = 25 nr = npts-2 np = npts+2 nt = npts N = 3*nr*np*nt num_steps = 2000 !Number of smvps to do in loop value = 1.0 IF (rank .eq. 0) THEN print*,'---------------------------------------------------' print*,'Running MAS PETSc tests with:' print*,'nr: ',nr print*,'np: ',np print*,'nt: ',nt print*,'N: ',N print*,'num_steps: ',num_steps print*,'MPI cores: ',size print*,'---------------------------------------------------' END IF print*,'MPI TEST: My rank is:',rank !Create x, y, and b vectors: call VecCreate(PETSC_COMM_WORLD,x,ierr) call VecSetType(x, VECMPI,ierr) call VecSetSizes(x,PETSC_DECIDE,N,ierr) call VecSetFromOptions(x,ierr) call VecCreate(PETSC_COMM_WORLD,y,ierr) call VecSetType(y, VECMPI,ierr) call VecSetSizes(y,PETSC_DECIDE,N,ierr) call VecSetFromOptions(y,ierr) call VecCreate(PETSC_COMM_WORLD,b,ierr) call VecSetType(b, VECMPI,ierr) call VecSetSizes(b,PETSC_DECIDE,N,ierr) call VecSetFromOptions(b,ierr) !Initialize y vector to 0: value = 0.0 call VecSet(y,value,ierr) !Initialize b and x vectors to (different) random values: IF (rank .eq. 0) THEN allocate(xvec(N)) allocate(yvec(N)) allocate(veci(N)) !Create list of indices (0-based for petsc) DO i=0,N-1 veci(i+1) = i END DO call RANDOM_NUMBER(xvec) call VecSetValues(b,N,veci,xvec,INSERT_VALUES,ierr) call RANDOM_NUMBER(xvec) call VecSetValues(x,N,veci,xvec,INSERT_VALUES,ierr) deallocate(veci) END IF !Assemble the x and b vectors across all cores: call VecAssemblyBegin(x,ierr); call VecAssemblyEnd(x,ierr); call VecAssemblyBegin(y,ierr); call VecAssemblyEnd(y,ierr); call VecAssemblyBegin(b,ierr); call VecAssemblyEnd(b,ierr); !call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr) !call VecView(y,PETSC_VIEWER_STDOUT_WORLD,ierr) !Create matrix: call MatCreate(PETSC_COMM_WORLD,A,ierr) call MatSetType(A,MATMPIAIJ,ierr) call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N,ierr) call MatSetFromOptions(A,ierr) !print*,'3nrt: ',3*nr*nt i = 16 IF(size .eq. 1) THEN j = 0 ELSE j = 8 END IF call MatMPIAIJSetPreallocation(A,i,PETSC_NULL_INTEGER, & j,PETSC_NULL_INTEGER,ierr) !Do not call this if using preallocation! !call MatSetUp(A,ierr) call MatGetOwnershipRange(A,istart,iend,ierr) print*,'Rank ',rank,' has rows ',istart,' to ',iend !Get MAS matrix in CSR format (random numbers for now): call GET_RAND_MAS_MATRIX(CSR_A,CSR_AI,CSR_AJ,nr,nt,np,M) IF (rank .eq. 0) THEN print*,'Number of non-zero entries in matrix: ',M !print*,'sss',CSR_A(:) !call csr Ax for validation/timing write(*,'(a)',advance='yes') & ' Computing y=Ax with RON_CSR_AX...' ttmp1 = etime(elapsed) DO i=1,num_steps CALL RON_CSR_AX(CSR_A,CSR_AI,CSR_AJ,xvec,yvec,N,M) END DO ttmp2 = etime(elapsed) CSRtime = 3.e9*(ttmp2-ttmp1)/(1.0*N*num_steps) !CSRMFLOPS = ((2.0*M*num_steps)/(ttmp2-ttmp1))/1024/1024/1024 write(*,'(a)',advance='yes') ' ...done!' print*,'||y||: ',norm(yvec,N) print*,'y(5)= ',yvec(5) !Store matrix values one-by-one (inefficient: better way ! more complicated - implement later) print*,'Storing MAS matrix into PETSc matrix...' END IF DO i=istart+1,iend !print*,'numofnonzerosinrowi:',CSR_AJ(i+1)-CSR_AJ(i)+1 DO j=CSR_AJ(i)+1,CSR_AJ(i+1) call MatSetValue(A,i-1,CSR_AI(j),CSR_A(j), & INSERT_VALUES,ierr) END DO END DO IF(rank .eq. 0) THEN print*,'...done!' END IF !Assemble matrix A across all cores: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) print*,'rank ',rank,' about to call MatAssemblyEnd()...' call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) !call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr) !call MatView(A,PETSC_VIEWER_DRAW_WORLD,ierr) IF (rank .eq. 0) THEN !Clean up CSR matrices: deallocate(CSR_A) deallocate(CSR_AI) deallocate(CSR_AJ) END IF !Call matrix-vector-product y=Ax num-steps times: IF (rank .eq. 0) THEN write(*,'(a)',advance='yes') & ' Computing y=Ax with PETSc MatMult...' ttmp1 = etime(elapsed) END IF DO i=1,num_steps call MatMult(A,x,y,ierr) END DO IF (rank .eq. 0) THEN ttmp2 = etime(elapsed) flops = ((2.0*M*num_steps)/(ttmp2-ttmp1))/1024/1024/1024 PETSctime = 3.e9*(ttmp2-ttmp1)/(1.0*N*num_steps) END IF !Get 2-norm of result vector using petsc: call VecNorm(y,NORM_2,value,ierr) IF (rank .eq. 0) THEN write(*,'(a)',advance='yes') ' ...done!' print*,'||y||: ',value call VecGetValues(y,1,4,value,ierr) print*,'y(5)= ',value print*,'RON_CSR_AX y=Ax time: ',CSRtime,'nsec/mp.' print*,'PETSc y=Ax time: ',PETSctime,'nsec/mp.' print*,'PETSc y=Ax flops: ',flops,'GFLOPS.' END IF !call VecView(y,PETSC_VIEWER_STDOUT_WORLD,ierr) c--------------------KSP LINEAR SOLVER TEST----------------------- !Initialize solver call KSPCreate(PETSC_COMM_WORLD,cgsolve,ierr) call KSPSetOperators(cgsolve,A,A,SAME_NONZERO_PATTERN,ierr) call KSPSetFromOptions(cgsolve,ierr) !Call CG with preconditioner set by input flag: !call KSPSolve(cgsolve,b,x,ierr); IF(rank .eq. 0) THEN !Get number of iterations: call KSPGetIterationNumber(cgsolve,i,ierr); call KSPGetConvergedReason(cgsolve,reason) IF (reason > 0) THEN print*,'Converged in ',i,' iterations' call KSPGetSolution(cgsolve,x,ierr); ELSE print*,'Did not converge! Number of iteratioins: ',i END IF END IF !Must call this again before re-solving if A changes: !call KSPSetOperators(cgsolve,A,A,SAME_NONZERO_PATTERN,ierr) call KSPDestroy(cgsolve,ierr) !Set-up Shell Matrix: !Clean up: IF (rank .eq. 0) THEN deallocate(xvec) deallocate(yvec) END IF call MatDestroy(A,ierr) call VecDestroy(x,ierr) call VecDestroy(y,ierr) call VecDestroy(b,ierr) call PetscFinalize(ierr) c----------------------------------------------------------------- c----------------------------------------------------------------- c----------------------------------------------------------------- contains !Generate MAS sparsity pattern of random numbers in DIA format: SUBROUTINE GET_RAND_MAS_MATRIX(CSR_A,CSR_AI,CSR_AJ,nr,nt,np,M) implicit none real*8, dimension(:,:),allocatable :: DIA_A integer, dimension(:), allocatable :: DIA_AI real*8, dimension(:), allocatable :: CSR_A integer, dimension(:), allocatable :: CSR_AI integer, dimension(:), allocatable :: CSR_AJ integer :: N,M,ndiags,ai,i,j,k,istart,iend,nr,nt,np N = 3*nr*np*nt ndiags = 29 allocate(DIA_AI(ndiags)) !Set positions of diags: DIA_AI(1) = -3*nr*nt DIA_AI(2) = -3*nr*nt + 1 DIA_AI(3) = -3*nr*nt + 2 DIA_AI(4) = -3*(nr*nt-1) + 2 DIA_AI(5) = -3*(nr*nt-nr) + 1 DIA_AI(6) = -3*nr - 1 DIA_AI(7) = -3*nr DIA_AI(8) = -3*nr + 1 DIA_AI(9) = -3*(nr-1) + 1 DIA_AI(10) = -5 DIA_AI(11) = -4 DIA_AI(12) = -3 DIA_AI(13) = -2 DIA_AI(14) = -1 DIA_AI(15) = 0 !Assume symmetric diag-structure (not values): DO i=(ndiags-ndiags/2),ndiags DIA_AI(i) = -DIA_AI(ndiags+1-i) END DO allocate(DIA_A(N,ndiags)) DIA_A(:,:) = -9999 CALL RANDOM_SEED() !MAS sparse pattern DO k=1,ndiags IF(DIA_AI(k) < 0) THEN istart = -DIA_AI(k)+1 iend = N ELSE istart = 1 iend = N-DIA_AI(k) END IF !Seven full diagonals: IF(k .eq. 15) THEN call RANDOM_NUMBER(DIA_A(istart:iend,k)) END IF IF((k .eq. 12) .or. (k .eq. 18)) THEN call RANDOM_NUMBER(DIA_A(istart:iend,k)) END IF IF((k .eq. 7) .or. (k .eq. 23)) THEN call RANDOM_NUMBER(DIA_A(istart:iend,k)) END IF IF((k .eq. 1) .or. (k .eq. 29)) THEN call RANDOM_NUMBER(DIA_A(istart:iend,k)) END IF !Two fill-2, skip 1 IF((k .eq. 14) .or. (k .eq. 16)) THEN i = istart DO CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT END DO END IF !Six fill1-skip2 IF((k .eq. 13) .or. (k .eq. 17)) THEN i = istart DO CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT END DO END IF IF((k .eq. 11) .or. (k .eq. 19)) THEN i = istart DO CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT END DO END IF IF((k .eq. 10) .or. (k .eq. 20)) THEN i = istart DO CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT END DO END IF !Six skipinit-fill1-skip2 IF((k .eq. 9) .or. (k .eq. 21)) THEN i = istart DIA_A(i,k) = 0 i = i+1 DO CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT END DO END IF IF((k .eq. 8) .or. (k .eq. 22)) THEN i = istart DIA_A(i,k) = 0 i = i+1 DO CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT END DO END IF IF((k .eq. 6) .or. (k .eq. 24)) THEN i = istart DIA_A(i,k) = 0 i = i+1 DO CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT END DO END IF !Eight skip2-fill1 IF((k .eq. 5) .or. (k .eq. 25)) THEN i = istart DO DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT END DO END IF IF((k .eq. 4) .or. (k .eq. 26)) THEN i = istart DO DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT END DO END IF IF((k .eq. 3) .or. (k .eq. 27)) THEN i = istart DO DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT END DO END IF IF((k .eq. 2) .or. (k .eq. 28)) THEN i = istart DO DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT DIA_A(i,k) = 0 i = i+1 IF(i .gt. iend) EXIT CALL RANDOM_NUMBER(DIA_A(i,k)) i = i+1 IF(i .gt. iend) EXIT END DO END IF END DO M = 0 DO j=1,ndiags DO i=1,N IF((DIA_A(i,j) .ne. -9999) .and. (DIA_A(i,j) .ne. 0)) THEN M = M+1 END IF END DO END DO allocate(CSR_A(M)) allocate(CSR_AI(M)) allocate(CSR_AJ(N+1)) !Use zero-indexing for A for use with petsc: ai=1 CSR_AJ(1) = 0 DO i=1,N DO k=1,ndiags IF((DIA_A(i,k) .ne. 0) .and. (DIA_A(i,k) .ne. -9999)) THEN CSR_AI(ai) = i + DIA_AI(k) - 1 CSR_A(ai) = DIA_A(i,k) ai = ai+1 END IF END DO CSR_AJ(i+1) = ai-1 END DO deallocate(DIA_A,DIA_AI) END SUBROUTINE c------------------------------------------------------------ !Compute y=Ax with CSR matrix using zero-indexing in A SUBROUTINE RON_CSR_AX(CSR_A,CSR_AI,CSR_AJ,xvec,yvec,N,M) real*8, dimension(M) :: CSR_A integer, dimension(M) :: CSR_AI integer, dimension(N+1) :: CSR_AJ real*8, dimension(N) :: xvec real*8, dimension(N) :: yvec integer :: N,M,i,j real*8 :: sum DO i=1,N sum = 0.0 DO j=CSR_AJ(i)+1, CSR_AJ(i+1) sum = sum + CSR_A(j)*xvec(CSR_AI(j)+1) END DO yvec(i) = sum END DO END SUBROUTINE c------------------------------------------------------------ !Get 2-norm of matrix x: FUNCTION norm (x,N) implicit none real*8, dimension(N) :: x real*8 :: norm integer :: i,N norm = 0.0 DO i=1,N norm = norm + x(i)**2 END DO norm = SQRT(norm) RETURN END FUNCTION END PROGRAM