The "PETSC_COMM_SELF" worked like a charm for me, thanks alot
On Mon, Mar 10, 2014 at 4:50 PM, Matthew Knepley <[email protected]> wrote: > On Mon, Mar 10, 2014 at 9:36 AM, Mohammad Bahaa > <[email protected]>wrote: > >> Hi, >> >> I'm pretty new to PETSC, so pardon me if the question is primitive >> somehow, I used *METIS *to partition my grid (represented by a system of >> linear equations Ax=b) to a number of sub-grids, say 4 sub-grids, with 4 >> different systems of linear Equations (A1x1=b1, A2x2=b2, ...), can anyone >> post an example showing how to solve these "n" sub-systems simultaneously, >> I've tried the following program, but it's not working correctly, as when I >> use *MatGetOwnershipRange *in each process I find that A1 ownership >> range is 1/4 the matrix size for the first process, while it should be all >> of it. >> > > I will answer this two ways. First, here is the "PETSc strategy" for doing > the same thing. > > 1) Write a code that assembles and solves the entire thing > > 2) Use -pc_type bjacobi -ksp_type preonly -sub_ksp_type <whatever you > want> -mat > > 3) You can use ParMetis inside PETSc with the MatPartitioning > > This will solve the individual systems with no coupling (this is what it > sounds like you want above). > > If you want to manage everything yourself, and you want to form individual > systems on every process, > just create the solvers using a smaller communicator. PETSC_COMM_SELF > means that every system > is serial. You can make smaller comms with MPI_Comm_split() if you want > smaller comms, but some > parallelism for each system. > > Thanks, > > Matt > > >> subroutine test_drive_2 >> >> implicit none >> >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! Include files >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! >> ! This program uses CPP for preprocessing, as indicated by the use of >> ! PETSc include files in the directory petsc/include/finclude. This >> ! convention enables use of the CPP preprocessor, which allows the use >> ! of the #include statements that define PETSc objects and variables. >> ! >> ! Use of the conventional Fortran include statements is also supported >> ! In this case, the PETsc include files are located in the directory >> ! petsc/include/foldinclude. >> ! >> ! Since one must be very careful to include each file no more than once >> ! in a Fortran routine, application programmers must exlicitly list >> ! each file needed for the various PETSc components within their >> ! program (unlike the C/C++ interface). >> ! >> ! See the Fortran section of the PETSc users manual for details. >> ! >> ! The following include statements are required for KSP Fortran programs: >> ! petscsys.h - base PETSc routines >> ! petscvec.h - vectors >> ! petscmat.h - matrices >> ! petscksp.h - Krylov subspace methods >> ! petscpc.h - preconditioners >> ! Other include statements may be needed if using additional PETSc >> ! routines in a Fortran program, e.g., >> ! petscviewer.h - viewers >> ! petscis.h - index sets >> ! >> >> ! main includes >> #include <finclude/petscsys.h> >> #include <finclude/petscvec.h> >> #include <finclude/petscmat.h> >> #include <finclude/petscksp.h> >> #include <finclude/petscpc.h> >> >> ! includes for F90 specific functions >> #include <finclude/petscvec.h90> >> #include <finclude/petscmat.h90> >> #include <finclude/petscksp.h90> >> #include <finclude/petscpc.h90> >> ! >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! Variable declarations >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! >> ! Variables: >> ! ksp - linear solver context >> ! ksp - Krylov subspace method context >> ! pc - preconditioner context >> ! x, b, u - approx solution, right-hand-side, exact solution vectors >> ! A - matrix that defines linear system >> ! its - iterations for convergence >> ! norm - norm of error in solution >> ! >> Vec x,b,u >> Mat A >> KSP ksp >> PC pc >> PetscReal norm,tol >> PetscErrorCode ierr >> PetscInt i,n,col(3),its,i1,i2,i3 >> PetscBool flg >> PetscMPIInt size,rank >> PetscScalar none,one,value(3), testa >> >> PetscScalar, pointer :: xx_v(:) >> PetscScalar, allocatable, dimension(:) :: myx >> PetscOffset i_x >> >> !real(4) :: myx(10), myu(10), myb(10) >> !real(8), allocatable, dimension(:) :: myx >> >> integer :: ic, nc, ncmax, nz, ncols, j >> integer :: fileunit, ione >> integer, allocatable, dimension(:,:) :: neighb, cols >> integer, allocatable, dimension(:) :: nnz, vcols >> real(8), allocatable, dimension(:,:) :: acoef, vals >> real(8), allocatable, dimension(:) :: ap, su, vvals >> character (len=100) :: rankstring, filename, folder >> >> real(8) :: atol, rtol, dtol >> integer :: mxit, istart, iend >> real(8) :: rvar, minvalx >> >> >> call PetscInitialize(PETSC_NULL_CHARACTER,ierr) >> call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr) >> call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) >> >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! Load the linear system >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> >> !mbs read data file for experimentation >> !if(rank.EQ.0) read(*,'(A)'), folder >> folder = '60' >> write(rankstring,'(I)'), rank >> filename = trim(adjustl(folder)) // '/linearsys_' // >> trim(adjustl(rankstring)) // '.txt' >> fileunit = 9000 + rank >> open(unit=fileunit, file=trim(filename)) >> >> read(fileunit,*), ncmax >> read(fileunit,*), nc >> read(fileunit,*), >> >> allocate( neighb(6,ncmax), acoef(6,ncmax), ap(ncmax), su(ncmax) ) >> allocate( nnz(0:ncmax-1), cols(6,0:ncmax-1), vals(6,0:ncmax-1) ) >> allocate( vcols(0:ncmax-1), vvals(0:ncmax-1) ) >> !allocate( xx_v(ncmax), myx(ncmax) ) >> !allocate( xx_v(0:ncmax), myx(0:ncmax) ) >> allocate( myx(ncmax) ) >> >> do i=1,nc >> read(fileunit,'(I)'), ic >> read(fileunit,'(6I)'), ( neighb(j,ic), j=1,6 ) >> read(fileunit,'(6F)'), ( acoef(j,ic), j=1,6 ) >> read(fileunit,'(2F)'), ap(ic), su(ic) >> read(fileunit,*), >> enddo >> >> close(fileunit) >> >> nz = 7 >> nnz = 0 >> do ic=0,nc-1 >> >> ! values for coefficient matrix (diagonal) >> nnz(ic) = nnz(ic) + 1 >> cols(nnz(ic),ic) = ic >> vals(nnz(ic),ic) = ap(ic+1) >> >> ! values for coefficient matrix (off diagonal) >> do j=1,6 >> if(neighb(j,ic+1).GT.0)then >> nnz(ic) = nnz(ic) + 1 >> cols(nnz(ic),ic) = neighb(j,ic+1) - 1 >> vals(nnz(ic),ic) = acoef(j,ic+1) >> endif >> enddo >> >> ! values for RHS >> vcols(ic) = ic >> vvals(ic) = su(ic+1) >> >> enddo >> >> ! add dummy values for remaining rows (if any) >> if(ncmax.GT.nc)then >> do ic=nc,ncmax-1 >> ! coeff matrix >> nnz(ic) = 1 >> cols(nnz(ic),ic) = ic >> vals(nnz(ic),ic) = 1.0 >> ! RHS >> vcols(ic) = ic >> vvals(ic) = 0.0 >> enddo >> endif >> >> print*, 'rank', rank, 'says nc is', nc >> >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! Beginning of program >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> >> !if (size .ne. 1) then >> ! call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) >> ! if (rank .eq. 0) then >> ! write(6,*) 'This is a uniprocessor example only!' >> ! endif >> ! SETERRQ(PETSC_COMM_WORLD,1,' ',ierr) >> !endif >> >> ione = 1 >> none = -1.0 >> one = 1.0 >> n = ncmax >> i1 = 1 >> i2 = 2 >> i3 = 3 >> call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr) >> >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! Compute the matrix and right-hand-side vector that define >> ! the linear system, Ax = b. >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> >> ! Create matrix. When using MatCreate(), the matrix format can >> ! be specified at runtime. >> >> call MatCreate(PETSC_COMM_WORLD,A,ierr) >> !call MatCreateSeqAij(PETSC_COMM_SELF,A,ierr) >> call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr) >> call MatSetFromOptions(A,ierr) >> call MatSetUp(A,ierr) >> >> call MatGetOwnershipRange(A,istart,iend,ierr) >> print*, rank, istart, iend >> >> ! Assemble matrix. >> ! - Note that MatSetValues() uses 0-based row and column numbers >> ! in Fortran as well as in C (as set here in the array "col"). >> >> ! value(1) = -1.0 >> ! value(2) = 2.0 >> ! value(3) = -1.0 >> ! do 50 i=1,n-2 >> ! col(1) = i-1 >> ! col(2) = i >> ! col(3) = i+1 >> ! call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr) >> !50 continue >> ! i = n - 1 >> ! col(1) = n - 2 >> ! col(2) = n - 1 >> ! call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr) >> ! i = 0 >> ! col(1) = 0 >> ! col(2) = 1 >> ! value(1) = 2.0 >> ! value(2) = -1.0 >> ! call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr) >> ! call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) >> ! call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) >> >> do ic=0,ncmax-1 >> call >> MatSetValues(A,ione,ic,nnz(ic),cols(1:nnz(ic),ic),vals(1:nnz(ic),ic),INSERT_VALUES,ierr) >> enddo >> >> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) >> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) >> >> ! Create vectors. Note that we form 1 vector from scratch and >> ! then duplicate as needed. >> >> call VecCreate(PETSC_COMM_WORLD,x,ierr) >> !call VecCreateSeq(PETSC_COMM_SELF,x,ierr) >> call VecSetSizes(x,PETSC_DECIDE,n,ierr) >> call VecSetFromOptions(x,ierr) >> call VecDuplicate(x,b,ierr) >> call VecDuplicate(x,u,ierr) >> >> ! Set exact solution; then compute right-hand-side vector. >> >> call VecSet(u,one,ierr) >> call VecSet(x,one*10,ierr) >> !call MatMult(A,u,b,ierr) >> >> ! set source terms vector >> call VecSetValues(b,ncmax,vcols,vvals,INSERT_VALUES,ierr) >> call VecAssemblyBegin(b,ierr) >> call VecAssemblyEnd(b,ierr) >> >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! Create the linear solver and set various options >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> >> ! Create linear solver context >> >> call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) >> >> ! Set operators. Here the matrix that defines the linear system >> ! also serves as the preconditioning matrix. >> >> call KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN,ierr) >> >> ! Set linear solver defaults for this problem (optional). >> ! - By extracting the KSP and PC contexts from the KSP context, >> ! we can then directly directly call any KSP and PC routines >> ! to set various options. >> ! - The following four statements are optional; all of these >> ! parameters could alternatively be specified at runtime via >> ! KSPSetFromOptions(); >> >> call KSPGetPC(ksp,pc,ierr) >> call PCSetType(pc,PCJACOBI,ierr) >> atol = 1.d-12 >> rtol = 1.d-12 >> dtol = 1.d10 >> mxit = 100 >> ! call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION, & >> ! & PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr) >> >> call KSPSetTolerances(ksp,atol,rtol,dtol,mxit,ierr) >> >> ! Set runtime options, e.g., >> ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> >> ! These options will override those specified above as long as >> ! KSPSetFromOptions() is called _after_ any other customization >> ! routines. >> >> call KSPSetFromOptions(ksp,ierr) >> >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! Solve the linear system >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> >> call KSPSetType(ksp,KSPBCGS,ierr) >> >> call KSPSolve(ksp,b,x,ierr) >> >> ! View solver info; we could instead use the option -ksp_view >> >> !call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr) >> >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> ! Check solution and clean up >> ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> >> call VecGetArrayF90(x,xx_v,ierr) >> !xx_v = 5.1d0 >> !do ic=1,ncmax >> ! myx(ic) = xx_v(ic) >> !enddo >> myx = xx_v >> !call VecGetArray(x,myx,i_x,ierr) >> !value = x_array(i_x + 1) >> !call VecRestoreArray(x,myx,i_x,ierr) >> !rvar = xx_v(3) >> call VecRestoreArrayF90(x,xx_v,ierr) >> ! >> !call VecGetArrayF90(b,xx_v,ierr) >> !myb = xx_v >> !call VecRestoreArrayF90(x,xx_v,ierr) >> ! >> !call VecView(x,PETSC_VIEWER_STDOUT_SELF) >> !call MatView(a,PETSC_VIEWER_STDOUT_SELF) >> >> !print*, 'rank', rank, 'says max x is', maxval(myx) >> ! print*, xx_v >> >> ! Check the error >> >> call MatMult(A,x,u,ierr) >> call VecAXPY(u,none,b,ierr) >> call VecNorm(u,NORM_2,norm,ierr) >> call KSPGetIterationNumber(ksp,its,ierr) >> >> if (norm .gt. 1.e-12) then >> write(6,100) norm,its >> else >> write(6,200) its >> endif >> 100 format('Norm of error = ',e11.4,', Iterations = ',i5) >> 200 format('Norm of error < 1.e-12,Iterations = ',i5) >> >> !call KSPGetSolution(ksp,myx) >> >> minvalx = 1.0e15 >> do ic=1,ncmax >> if(myx(ic).LT.minvalx) minvalx = myx(ic) >> enddo >> >> !write(*,300), rank, maxval(myx(1:nc)), minvalx >> >> if(rank.EQ.0) print*, myx >> >> 300 format('Rank ', I1, ' says max/min x are: ', F12.4, ' / ', F12.4) >> >> ! Free work space. All PETSc objects should be destroyed when they >> ! are no longer needed. >> >> call VecDestroy(x,ierr) >> call VecDestroy(u,ierr) >> call VecDestroy(b,ierr) >> call MatDestroy(A,ierr) >> call KSPDestroy(ksp,ierr) >> call PetscFinalize(ierr) >> >> continue >> >> end >> >> -- >> Mohamamd Bahaa ElDin >> > > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener > -- Mohamamd Bahaa ElDin
