! ! Modified from ex15f.F for testing Superlu ! Solves a linear system in parallel with KSP. ! Current only work with 1 processors ! ------------------------------------------------------------------------- program main implicit none ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Include files ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #include #include #include #include #include #include #include ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Variable declarations ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Vec x,b,u Mat A Mat F PC pc KSP ksp PetscScalar v,one,neg_one,zero PetscReal norm,tol PetscErrorCode ierr PetscInt i,j,II,JJ,Istart PetscInt Iend,m,n,l,i1,its,five PetscInt mm,isol PetscMPIInt rank PetscBool flg,flg_ksp,initialguess PetscViewer fd Character(len=256) filename_a Character(len=256) filename_b Character(len=256) filename_c Character(len=256) filename_d Character(len=256) foldermatrix flg = PETSC_FALSE initialguess = PETSC_FALSE ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Beginning of program ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call PetscInitialize(PETSC_NULL_CHARACTER,ierr) one = 1.0 neg_one = -1.0 zero = 0.0 call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) ! Determine files matrix and right-hand-side vector call PetscOptionsGetString(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, & "-f0",filename_a,flg,ierr) CHKERRQ(ierr) if (flg) then ! Open binary file. Note that we use FILE_MODE_READ to indicate ! reading from this file. if (rank .eq. 0) then write(*,*) "-->loac matrix a" end if call PetscViewerBinaryOpen(PETSC_COMM_WORLD, & trim(filename_a), FILE_MODE_READ,fd,ierr) CHKERRQ(ierr) ! Load the matrix and vector; then destroy the viewer. call MatCreate(PETSC_COMM_WORLD,A,ierr) CHKERRQ(ierr) call MatSetFromOptions(A,ierr) CHKERRQ(ierr) call MatLoad(A,fd,ierr) CHKERRQ(ierr) call MatSetFromOptions(A,ierr) call PetscOptionsInsertString(PETSC_NULL_OBJECT, & "-mat_superlu_dist_fact SamePattern",ierr) endif ! Load right-hand-side vector flg = PETSC_FALSE call PetscOptionsGetString(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, & "-rhs",filename_b,flg,ierr) CHKERRQ(ierr) call VecCreate(PETSC_COMM_WORLD,b,ierr) CHKERRQ(ierr) if (flg) then if (rank .eq. 0) then write(*,*) "-->load rhs b" end if call PetscViewerDestroy(fd,ierr) CHKERRQ(ierr) call PetscViewerBinaryOpen(PETSC_COMM_WORLD, & trim(filename_b),FILE_MODE_READ,fd,ierr) CHKERRQ(ierr) call VecSetFromOptions(b,ierr) CHKERRQ(ierr) call VecLoad(b,fd,ierr) CHKERRQ(ierr) endif call PetscViewerDestroy(fd,ierr) ! Initialize solution vector x and u call MatCreateVecs(A,x,PETSC_NULL_OBJECT,ierr) CHKERRQ(ierr) call VecDuplicate(b,u,ierr) CHKERRQ(ierr) ! Check size of matrix and vector call MatGetLocalSize(A,m,n,ierr) CHKERRQ(ierr) call MatGetSize(A,mm,PETSC_NULL_INTEGER,ierr) CHKERRQ(ierr) call VecGetSize(x,l,ierr) CHKERRQ(ierr) if (rank .eq. 0) then write(*,*) "size l,m,n,mm",l,m,n,mm end if ! Set solution to zero call VecSet(x,zero,ierr) CHKERRQ(ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Create the linear solver and set various options ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) CHKERRQ(ierr) call KSPSetInitialGuessNonzero(ksp,initialguess,ierr) CHKERRQ(ierr) call KSPSetOperators(ksp,A,A,ierr) CHKERRQ(ierr) ! Set tolerance tol = 1.e-5 call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL, & PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr) CHKERRQ(ierr) ! Test SUPERLU_DIST flg_ksp = PETSC_FALSE call PetscOptionsGetBool(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, & "-superlu_default",flg_ksp,PETSC_NULL_CHARACTER,ierr) CHKERRQ(ierr) if (flg_ksp) then call KSPSetType(ksp,KSPPREONLY,ierr) CHKERRQ(ierr) call KSPGetPC(ksp,pc,ierr) CHKERRQ(ierr) call PCSetType(pc,PCLU,ierr) CHKERRQ(ierr) call PCFactorSetMatSolverPackage(pc,MATSOLVERSUPERLU_DIST,ierr) CHKERRQ(ierr) call PCFactorSetUpMatSolverPackage(pc,ierr) CHKERRQ(ierr) call PCFactorGetMatrix(pc,F,ierr) CHKERRQ(ierr) endif call KSPSetFromOptions(ksp,ierr) CHKERRQ(ierr) call KSPSetUp(ksp,ierr) CHKERRQ(ierr) call KSPSetUpOnBlocks(ksp,ierr) CHKERRQ(ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Solve the linear system ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !****************************************** ! First test, solver the first equation, the first equation is read from ! the input options !****************************************** call KSPSolve(ksp,b,x,ierr) CHKERRQ(ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Check solution and clean up ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call MatMult(A,x,u,ierr) CHKERRQ(ierr) call VecAXPY(u,neg_one,b,ierr) CHKERRQ(ierr) call VecNorm(u,NORM_2,norm,ierr) CHKERRQ(ierr) call KSPGetIterationNumber(ksp,its,ierr) CHKERRQ(ierr) if (rank .eq. 0) then if (norm .gt. 1.e-12) then write(6,100) norm,its else !write(6,110) its write(6,100) norm,its endif endif 100 format('Norm of error ',1pe11.4,' iterations ',i5) ! 110 format('Norm of error < 1.e-12,iterations ',i5) !****************************************** ! Second test, solver matrix 168 to 172, user need to specify ! the name and folder where these matrices and rhs are stored. !****************************************** flg = PETSC_FALSE call PetscOptionsGetString(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, & "-loop_matrices",filename_c,flg,ierr) flg = PETSC_FALSE call PetscOptionsGetString(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, & "-loop_folder",foldermatrix,flg,ierr) if (flg) then do isol = 168, 172 if (rank == 0) then write(*,*) "-->Test for matrix ",isol end if ! Load matrix write(filename_d, *) isol filename_d = trim(adjustl(foldermatrix))//"/a_"// & trim(filename_c)//"_"// & trim(adjustl(filename_d))//".bin" call MatZeroEntries(A, ierr) CHKERRQ(ierr) call PetscViewerBinaryOpen(PETSC_COMM_WORLD, & trim(filename_d), FILE_MODE_READ,fd,ierr) CHKERRQ(ierr) call MatSetFromOptions(A,ierr) CHKERRQ(ierr) call MatLoad(A,fd,ierr) CHKERRQ(ierr) call PetscViewerDestroy(fd,ierr) ! Load rhs write(filename_d, *) isol filename_d = "b_"//trim(filename_c)//"_"// & trim(adjustl(filename_d))//".bin" call VecSet(b,zero,ierr) CHKERRQ(ierr) call PetscViewerBinaryOpen(PETSC_COMM_WORLD, & trim(filename_b),FILE_MODE_READ,fd,ierr) CHKERRQ(ierr) call VecSetFromOptions(b,ierr) CHKERRQ(ierr) call VecLoad(b,fd,ierr) CHKERRQ(ierr) call PetscViewerDestroy(fd,ierr) ! Set x and u to zero call VecSet(x,zero,ierr) CHKERRQ(ierr) call VecSet(u,zero,ierr) CHKERRQ(ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Solve the linear system ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call KSPSolve(ksp,b,x,ierr) CHKERRQ(ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Check solution and clean up ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call MatMult(A,x,u,ierr) CHKERRQ(ierr) call VecAXPY(u,neg_one,b,ierr) CHKERRQ(ierr) call VecNorm(u,NORM_2,norm,ierr) CHKERRQ(ierr) call KSPGetIterationNumber(ksp,its,ierr) CHKERRQ(ierr) if (rank .eq. 0) then if (norm .gt. 1.e-12) then write(6,120) norm,its else !write(6,110) its write(6,120) norm,its endif endif 120 format('Norm of error ',1pe11.4,' iterations ',i5) end do end if ! Free work space. All PETSc objects should be destroyed when they ! are no longer needed. call KSPDestroy(ksp,ierr) CHKERRQ(ierr) call VecDestroy(u,ierr) CHKERRQ(ierr) call VecDestroy(x,ierr) CHKERRQ(ierr) call VecDestroy(b,ierr) CHKERRQ(ierr) call MatDestroy(A,ierr) CHKERRQ(ierr) if (rank .eq. 0) then write(*,*) "--> End of test, bye" end if ! Always call PetscFinalize() before exiting a program. call PetscFinalize(ierr) CHKERRQ(ierr) end