<div dir="ltr">The "<span style="font-family:arial,sans-serif;font-size:13px">PETSC_COMM_SELF" worked like a charm for me, thanks alot</span></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Mar 10, 2014 at 4:50 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div class="">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><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 class="h5"><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><font color="#888888"><div><br></div>-- <br><div dir="ltr">Mohamamd Bahaa ElDin</div>
</font></span></div></div>
</blockquote></div></div></div><span class="HOEnZb"><font color="#888888"><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
</font></span></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div dir="ltr">Mohamamd Bahaa ElDin</div>
</div>