<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Mon, Mar 10, 2014 at 9:36 AM, Mohammad Bahaa <span dir="ltr"><<a href="mailto:m.bahaa.eldein@gmail.com" target="_blank">m.bahaa.eldein@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">Hi,<div><br></div><div>I'm pretty new to PETSC, so pardon me if the question is primitive somehow, I used <b>METIS </b>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 <b>MatGetOwnershipRange </b>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.</div>
</div></blockquote><div><br></div><div>I will answer this two ways. First, here is the "PETSc strategy" for doing the same thing.</div><div><br></div><div>  1) Write a code that assembles and solves the entire thing</div>
<div><br></div><div>  2) Use -pc_type bjacobi -ksp_type preonly -sub_ksp_type <whatever you want> -mat</div><div><br></div><div>  3) You can use ParMetis inside PETSc with the MatPartitioning</div><div><br></div><div>
This will solve the individual systems with no coupling (this is what it sounds like you want above).</div><div><br></div><div>If you want to manage everything yourself, and you want to form individual systems on every process,</div>
<div>just create the solvers using a smaller communicator. PETSC_COMM_SELF means that every system</div><div>is serial. You can make smaller comms with MPI_Comm_split() if you want smaller comms, but some</div><div>parallelism for each system.</div>
<div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div dir="ltr"><div><div>subroutine test_drive_2</div><div>      </div><div>      implicit none</div><div><br></div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div>!                    Include files</div>

<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div>!</div><div>!  This program uses CPP for preprocessing, as indicated by the use of</div><div>!  PETSc include files in the directory petsc/include/finclude.  This</div>

<div>!  convention enables use of the CPP preprocessor, which allows the use</div><div>!  of the #include statements that define PETSc objects and variables.</div><div>!</div><div>!  Use of the conventional Fortran include statements is also supported</div>

<div>!  In this case, the PETsc include files are located in the directory</div><div>!  petsc/include/foldinclude.</div><div>!</div><div>!  Since one must be very careful to include each file no more than once</div><div>
!  in a Fortran routine, application programmers must exlicitly list</div>
<div>!  each file needed for the various PETSc components within their</div><div>!  program (unlike the C/C++ interface).</div><div>!</div><div>!  See the Fortran section of the PETSc users manual for details.</div><div>
!</div>
<div>!  The following include statements are required for KSP Fortran programs:</div><div>!     petscsys.h    - base PETSc routines</div><div>!     petscvec.h    - vectors</div><div>!     petscmat.h    - matrices</div><div>

!     petscksp.h    - Krylov subspace methods</div><div>!     petscpc.h     - preconditioners</div><div>!  Other include statements may be needed if using additional PETSc</div><div>!  routines in a Fortran program, e.g.,</div>

<div>!     petscviewer.h - viewers</div><div>!     petscis.h     - index sets</div><div>!</div><div><br></div><div>! main includes</div><div>#include <finclude/petscsys.h></div><div>#include <finclude/petscvec.h></div>

<div>#include <finclude/petscmat.h></div><div>#include <finclude/petscksp.h></div><div>#include <finclude/petscpc.h></div><div><br></div><div>! includes for F90 specific functions</div><div>#include <finclude/petscvec.h90></div>

<div>#include <finclude/petscmat.h90></div><div>#include <finclude/petscksp.h90></div><div>#include <finclude/petscpc.h90></div><div>!</div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>

<div>!                   Variable declarations</div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div>!</div><div>!  Variables:</div><div>!     ksp      - linear solver context</div>

<div>!     ksp      - Krylov subspace method context</div><div>!     pc       - preconditioner context</div><div>!     x, b, u  - approx solution, right-hand-side, exact solution vectors</div><div>!     A        - matrix that defines linear system</div>

<div>!     its      - iterations for convergence</div><div>!     norm     - norm of error in solution</div><div>!</div><div>      Vec              x,b,u</div><div>      Mat              A</div><div>      KSP              ksp</div>

<div>      PC               pc</div><div>      PetscReal        norm,tol</div><div>      PetscErrorCode   ierr</div><div>      PetscInt         i,n,col(3),its,i1,i2,i3</div><div>      PetscBool        flg</div><div>      PetscMPIInt      size,rank</div>

<div>      PetscScalar      none,one,value(3), testa</div><div>      </div><div>      PetscScalar, pointer :: xx_v(:)</div><div>      PetscScalar, allocatable, dimension(:) :: myx</div><div>      PetscOffset i_x</div><div>

      </div><div>      !real(4) :: myx(10), myu(10), myb(10)</div><div>      !real(8), allocatable, dimension(:) :: myx</div><div>      </div><div>      integer :: ic, nc, ncmax, nz, ncols, j</div><div>      integer :: fileunit, ione</div>

<div>      integer, allocatable, dimension(:,:) :: neighb, cols</div><div>      integer, allocatable, dimension(:)   :: nnz, vcols</div><div>      real(8), allocatable, dimension(:,:) :: acoef, vals</div><div>      real(8), allocatable, dimension(:)   :: ap, su, vvals</div>

<div>      character (len=100)                  :: rankstring, filename, folder</div><div>      </div><div>      real(8) :: atol, rtol, dtol</div><div>      integer :: mxit, istart, iend</div><div>      real(8) :: rvar, minvalx</div>

<div>      </div><div>      </div><div>      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)</div><div>      call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)</div><div>      call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)</div>

<div>      </div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div>!                 Load the linear system</div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>

<div>      </div><div>      !mbs read data file for experimentation</div><div>      !if(rank.EQ.0) read(*,'(A)'), folder</div><div>      folder = '60'</div><div>      write(rankstring,'(I)'), rank</div>

<div>      filename = trim(adjustl(folder)) // '/linearsys_' // trim(adjustl(rankstring)) // '.txt'</div><div>      fileunit = 9000 + rank</div><div>      open(unit=fileunit, file=trim(filename))</div><div>

      </div><div>      read(fileunit,*), ncmax</div><div>      read(fileunit,*), nc</div><div>      read(fileunit,*),</div><div>      </div><div>      allocate( neighb(6,ncmax), acoef(6,ncmax), ap(ncmax), su(ncmax) )</div>

<div>      allocate( nnz(0:ncmax-1), cols(6,0:ncmax-1), vals(6,0:ncmax-1) )</div><div>      allocate( vcols(0:ncmax-1), vvals(0:ncmax-1) )</div><div>      !allocate( xx_v(ncmax), myx(ncmax) )</div><div>      !allocate( xx_v(0:ncmax), myx(0:ncmax) )</div>

<div>      allocate( myx(ncmax) )</div><div>      </div><div>      do i=1,nc</div><div>          read(fileunit,'(I)'), ic</div><div>          read(fileunit,'(6I)'), ( neighb(j,ic), j=1,6 )</div><div>          read(fileunit,'(6F)'), ( acoef(j,ic), j=1,6 )</div>

<div>          read(fileunit,'(2F)'), ap(ic), su(ic)</div><div>          read(fileunit,*),</div><div>      enddo</div><div>      </div><div>      close(fileunit)</div><div>      </div><div>      nz = 7</div><div>
      nnz = 0</div>
<div>      do ic=0,nc-1</div><div>          </div><div>          ! values for coefficient matrix (diagonal)</div><div>          nnz(ic) = nnz(ic) + 1</div><div>          cols(nnz(ic),ic) = ic</div><div>          vals(nnz(ic),ic) = ap(ic+1)</div>

<div>          </div><div>          ! values for coefficient matrix (off diagonal)</div><div>          do j=1,6</div><div>              if(neighb(j,ic+1).GT.0)then</div><div>                  nnz(ic) = nnz(ic) + 1</div><div>

                  cols(nnz(ic),ic) = neighb(j,ic+1) - 1</div><div>                  vals(nnz(ic),ic) = acoef(j,ic+1)</div><div>              endif</div><div>          enddo</div><div>          </div><div>          ! values for RHS</div>

<div>          vcols(ic) = ic</div><div>          vvals(ic) = su(ic+1)</div><div>          </div><div>      enddo</div><div>      </div><div>      ! add dummy values for remaining rows (if any)</div><div>      if(<a href="http://ncmax.GT.nc" target="_blank">ncmax.GT.nc</a>)then</div>

<div>          do ic=nc,ncmax-1</div><div>              ! coeff matrix</div><div>              nnz(ic) = 1</div><div>              cols(nnz(ic),ic) = ic</div><div>              vals(nnz(ic),ic) = 1.0</div><div>              ! RHS</div>

<div>              vcols(ic) = ic</div><div>              vvals(ic) = 0.0</div><div>          enddo      </div><div>      endif</div><div>      </div><div>      print*, 'rank', rank, 'says nc is', nc</div>

<div><br></div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div>!                 Beginning of program</div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>

<div><br></div><div>      !if (size .ne. 1) then</div><div>      !   call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)</div><div>      !   if (rank .eq. 0) then</div><div>      !      write(6,*) 'This is a uniprocessor example only!'</div>

<div>      !   endif</div><div>      !   SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)</div><div>      !endif</div><div>      </div><div>      ione = 1</div><div>      none = -1.0</div><div>      one  = 1.0</div><div>      n    = ncmax</div>

<div>      i1 = 1</div><div>      i2 = 2</div><div>      i3 = 3</div><div>      call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr)</div><div><br></div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>

<div>!         Compute the matrix and right-hand-side vector that define</div><div>!         the linear system, Ax = b.</div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div><br></div>

<div>!  Create matrix.  When using MatCreate(), the matrix format can</div><div>!  be specified at runtime.</div><div><br></div><div>      call MatCreate(PETSC_COMM_WORLD,A,ierr)</div><div>      !call MatCreateSeqAij(PETSC_COMM_SELF,A,ierr)</div>

<div>      call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr)</div><div>      call MatSetFromOptions(A,ierr)</div><div>      call MatSetUp(A,ierr)</div><div>      </div><div>      call MatGetOwnershipRange(A,istart,iend,ierr)</div>

<div>      print*, rank, istart, iend</div><div><br></div><div>!  Assemble matrix.</div><div>!   - Note that MatSetValues() uses 0-based row and column numbers</div><div>!     in Fortran as well as in C (as set here in the array "col").</div>

<div><br></div><div>  !    value(1) = -1.0</div><div>  !    value(2) = 2.0</div><div>  !    value(3) = -1.0</div><div>  !    do 50 i=1,n-2</div><div>  !       col(1) = i-1</div><div>  !       col(2) = i</div><div>  !       col(3) = i+1</div>

<div>  !       call MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr)</div><div>  !50  continue</div><div>  !    i = n - 1</div><div>  !    col(1) = n - 2</div><div>  !    col(2) = n - 1</div><div>  !    call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)</div>

<div>  !    i = 0</div><div>  !    col(1) = 0</div><div>  !    col(2) = 1</div><div>  !    value(1) = 2.0</div><div>  !    value(2) = -1.0</div><div>  !    call MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr)</div><div>

  !    call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)</div><div>  !    call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)</div><div>  </div><div>       do ic=0,ncmax-1</div><div>           call MatSetValues(A,ione,ic,nnz(ic),cols(1:nnz(ic),ic),vals(1:nnz(ic),ic),INSERT_VALUES,ierr)</div>

<div>       enddo</div><div>       </div><div>       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)</div><div>       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)</div><div><br></div><div>!  Create vectors.  Note that we form 1 vector from scratch and</div>

<div>!  then duplicate as needed.</div><div><br></div><div>      call VecCreate(PETSC_COMM_WORLD,x,ierr)</div><div>      !call VecCreateSeq(PETSC_COMM_SELF,x,ierr)</div><div>      call VecSetSizes(x,PETSC_DECIDE,n,ierr)</div>

<div>      call VecSetFromOptions(x,ierr)</div><div>      call VecDuplicate(x,b,ierr)</div><div>      call VecDuplicate(x,u,ierr)</div><div><br></div><div>!  Set exact solution; then compute right-hand-side vector.</div>
<div>
<br></div><div>      call VecSet(u,one,ierr)</div><div>      call VecSet(x,one*10,ierr)</div><div>      !call MatMult(A,u,b,ierr)</div><div>      </div><div>      ! set source terms vector</div><div>      call VecSetValues(b,ncmax,vcols,vvals,INSERT_VALUES,ierr)</div>

<div>      call VecAssemblyBegin(b,ierr)</div><div>      call VecAssemblyEnd(b,ierr)</div><div>      </div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div>!          Create the linear solver and set various options</div>

<div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div><br></div><div>!  Create linear solver context</div><div><br></div><div>      call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)</div><div><br>

</div><div>!  Set operators. Here the matrix that defines the linear system</div><div>!  also serves as the preconditioning matrix.</div><div><br></div><div>      call KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN,ierr)</div>

<div><br></div><div>!  Set linear solver defaults for this problem (optional).</div><div>!   - By extracting the KSP and PC contexts from the KSP context,</div><div>!     we can then directly directly call any KSP and PC routines</div>

<div>!     to set various options.</div><div>!   - The following four statements are optional; all of these</div><div>!     parameters could alternatively be specified at runtime via</div><div>!     KSPSetFromOptions();</div>

<div><br></div><div>      call KSPGetPC(ksp,pc,ierr)</div><div>      call PCSetType(pc,PCJACOBI,ierr)</div><div>      atol = 1.d-12</div><div>      rtol = 1.d-12</div><div>      dtol = 1.d10</div><div>      mxit = 100</div>

<div>!      call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,     &</div><div>!     &     PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr)</div><div><br></div><div>       call KSPSetTolerances(ksp,atol,rtol,dtol,mxit,ierr)</div>

<div><br></div><div>!  Set runtime options, e.g.,</div><div>!      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol></div><div>!  These options will override those specified above as long as</div>

<div>!  KSPSetFromOptions() is called _after_ any other customization</div><div>!  routines.</div><div><br></div><div>      call KSPSetFromOptions(ksp,ierr)</div><div><br></div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>

<div>!                      Solve the linear system</div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div><br></div><div>      call KSPSetType(ksp,KSPBCGS,ierr)</div><div><br></div>

<div>      call KSPSolve(ksp,b,x,ierr)</div><div><br></div><div>!  View solver info; we could instead use the option -ksp_view</div><div><br></div><div>      !call KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD,ierr)</div><div><br>

</div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div>!                      Check solution and clean up</div><div>! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div>

<div><br></div><div>      call VecGetArrayF90(x,xx_v,ierr)</div><div>      !xx_v = 5.1d0</div><div>      !do ic=1,ncmax</div><div>      !    myx(ic) = xx_v(ic)</div><div>      !enddo</div><div>      myx = xx_v</div><div>
      !call VecGetArray(x,myx,i_x,ierr)</div>
<div>      !value = x_array(i_x + 1)</div><div>      !call VecRestoreArray(x,myx,i_x,ierr)</div><div>      !rvar = xx_v(3)</div><div>      call VecRestoreArrayF90(x,xx_v,ierr)</div><div>      !</div><div>      !call VecGetArrayF90(b,xx_v,ierr)</div>

<div>      !myb = xx_v</div><div>      !call VecRestoreArrayF90(x,xx_v,ierr)</div><div>      !</div><div>      !call VecView(x,PETSC_VIEWER_STDOUT_SELF)</div><div>      !call MatView(a,PETSC_VIEWER_STDOUT_SELF)</div><div>

      </div><div>      !print*, 'rank', rank, 'says max x is', maxval(myx)</div><div>     ! print*, xx_v</div><div>      </div><div>!  Check the error</div><div><br></div><div>      call MatMult(A,x,u,ierr)</div>

<div>      call VecAXPY(u,none,b,ierr)</div><div>      call VecNorm(u,NORM_2,norm,ierr)</div><div>      call KSPGetIterationNumber(ksp,its,ierr)</div><div>      </div><div>      if (norm .gt. 1.e-12) then</div><div>        write(6,100) norm,its</div>

<div>      else</div><div>        write(6,200) its</div><div>      endif</div><div> 100  format('Norm of error = ',e11.4,',  Iterations = ',i5)</div><div> 200  format('Norm of error < 1.e-12,Iterations = ',i5)</div>

<div>      </div><div>      !call KSPGetSolution(ksp,myx)</div><div>      </div><div>      minvalx = 1.0e15</div><div>      do ic=1,ncmax</div><div>          if(myx(ic).LT.minvalx) minvalx = myx(ic)</div><div>      enddo</div>

<div>      </div><div>      !write(*,300), rank, maxval(myx(1:nc)), minvalx</div><div>      </div><div>      if(rank.EQ.0) print*, myx</div><div>      </div><div>300 format('Rank ', I1, ' says max/min x are: ', F12.4, ' / ', F12.4)</div>

<div><br></div><div>!  Free work space.  All PETSc objects should be destroyed when they</div><div>!  are no longer needed.</div><div><br></div><div>      call VecDestroy(x,ierr)</div><div>      call VecDestroy(u,ierr)</div>

<div>      call VecDestroy(b,ierr)</div><div>      call MatDestroy(A,ierr)</div><div>      call KSPDestroy(ksp,ierr)</div><div>      call PetscFinalize(ierr)</div><div>      </div><div>      continue</div><div>      </div>

<div>end</div><span class=""><font color="#888888"><div><br></div>-- <br><div dir="ltr">Mohamamd Bahaa ElDin</div>
</font></span></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>