<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>