Thanks a lot Matthew. That was a great help.<div><div>It seems I was completely out of the track. I am looking at what you said, and made good progress til here...</div><div><br></div><div>Thanks again</div><div>Reza</div>
<div><br><br><div class="gmail_quote">On Wed, Dec 21, 2011 at 4:40 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com">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 class="im">On Wed, Dec 21, 2011 at 3:30 PM, Reza Madankan <span dir="ltr"><<a href="mailto:rm93@buffalo.edu" target="_blank">rm93@buffalo.edu</a>></span> wrote:<br></div><div class="gmail_quote"><div class="im">
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi Matthew;<div>Thank you for you quick reply.</div><div>Sorry for the wrong terminology, it's just few days I started working with Petsc. </div><div>The whole story is that I have a linear system Ax = b which I want to solve it using different preconditioners just to see the effect of each of those on the solution and residual error. I have already used PCJACOBI and PCILU to find the solution. I also would like to see the effect of using PCASM. PCASM on a single processor performs like PCILU, but its result will be different from PCILU if I implement it on multiple processors. So, I am trying to get it run for multiple processors case using the following command:<div>
<br>mpiexec -np 2 ./pcasm</div><div><br></div><div>but I get the error message that I showed before. Now my question is that how I can implement PCASM on multiple processors?</div></div></blockquote><div><br></div></div>
<div>It is not ASM that is the problem. You have a bug in your code. Use valgrind to find it.</div>
<div><br></div><div>The larger problem here is that you have completely misunderstood the point of PETSc.</div><div>Whatever you do, never change your code to use a different solver. You should play</div><div>with the PETSc examples for a while</div>
<div><br></div><div> cd src/snes/examples/tutorials</div><div> make ex5</div><div> ./ex5 -snes_monitor -ksp_monitor -snes_view</div><div><div> $MPIEXEC -n 2 ./ex5 -snes_monitor -ksp_monitor -snes_view</div></div><div>
<div><div> $MPIEXEC -n 2 ./ex5 -snes_monitor -ksp_monitor -snes_view -pc_type jacobi</div></div></div><div><div><div> $MPIEXEC -n 2 ./ex5 -snes_monitor -ksp_monitor -snes_view -pc_type asm</div></div></div><div><br></div>
<div> Matt</div><div><div class="h5"><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div>Thanks again</div><div>
Reza</div><div><br></div><div><br></div><div><br><div class="gmail_quote">On Wed, Dec 21, 2011 at 3:19 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>On Wed, Dec 21, 2011 at 2:12 PM, Reza Madankan <span dir="ltr"><<a href="mailto:rm93@buffalo.edu" target="_blank">rm93@buffalo.edu</a>></span> wrote:<br>
</div><div class="gmail_quote"><div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hello,<div>I am trying to parallelize the ASM preconditioner while solving a linear system. For more detail, I have printed relevant parts of code in the following:</div></blockquote><div><br></div></div><div>I have no idea what "parallelize ASM" might mean.</div>
<div><br></div><div>Below you likely have a bad value in b, the rhs you pass to the solver on line 487 of your code.</div><div><br></div><div> Matt</div><div><div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div><div><i>static char help[]="Reading in a matrix\n";</i></div>
<div><i><br></i></div><div><i> #include<stdio.h></i></div><div><i><br></i></div><div><i> #include<math.h></i></div><div><i><br></i></div><div><i> #include<stdlib.h></i></div><div><i><br></i></div><div><i> #include "petscvec.h" /* This enables us to use vectors. */</i></div>
<div><i><br></i></div><div><i> #include "petscmat.h" /* This enables us to use Matrices. It includes the</i></div><div><i> petscvec header file*/</i></div><div><i><br></i></div><div><i> #include "petscksp.h" /* Now we can solve linear systems. Solvers used are</i></div>
<div><i> KSP. */</i></div><div><i><br></i></div><div><i> extern PetscErrorCode MyKSPMonitor(KSP,PetscInt,PetscReal,void*);</i></div><div><i><br></i></div><div><i> int main(int argc, char **argv)</i></div><div><i><br></i></div>
<div><i> {</i></div><div><i><br></i></div><div><i> /* Declaration of Matrix A and some vectors x*/</i></div><div><i><br></i></div><div><i> Vec w, ypcq, mean_tmp, x, PzyTvec;</i></div><div><i> Mat Ym, Y, Mean, Ypcq, YpcqT, Pzz, InnProd, Para, Pzy, PzyT,PzzT,K, KT, PK, PzzTMP, Pmat, XYmat, KY;</i></div>
<div><i><br></i></div><div><i> FILE *fp, *fp1, *fp2, *fp3;</i></div><div><i><br></i></div><div><i> PetscInt index,i,j,k, ns=16, nw=100, tindex_f=120, col,tind, mz, nz, its;</i></div><div><i>/* </i></div><div><i>printf("Enter number of sensors (ns):\n");</i></div>
<div><i>scanf("%d",&ns);</i></div><div><i><br></i></div><div><i>printf("Enter number of time steps (tindex_f):\n");</i></div><div><i>scanf("%d",&tindex_f);</i></div><div><i>*/</i></div>
<div><i> PetscMPIInt<span style="white-space:pre-wrap"> </span>size;</i></div><div><i><br></i></div><div><i> PC pc;</i></div><div><i> KSP ksp;</i></div><div><i> PetscReal<span style="white-space:pre-wrap"> </span>norm, tol=1.e-14;</i></div>
<div><i> PetscBool<span style="white-space:pre-wrap"> </span>nonzeroguess = PETSC_FALSE;</i></div><div><i> PetscViewer<span style="white-space:pre-wrap"> </span>mviewer;</i></div><div><i><br>
</i></div><div><i><br></i></div><div><i> PetscScalar scalar,rhsy[ns*tindex_f*nw],rhsym[ns*tindex_f], rhsw[nw], Ymat[tindex_f*ns][nw], tmp, rhsp[2*nw], xvec[ns*tindex_f], pmat[4], xymat[2], rnorm;</i></div><div><i><br></i></div>
<div><i> PetscErrorCode ierr;</i></div><div><i><br></i></div><div><i> /* This part is needed for the help flag supplied at run-time*/</i></div><div><i><br></i></div><div><i> ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);</i></div>
<div><i> </i></div><div><i> ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);</i></div><div><i> //if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example");</i></div><div><i>PetscPrintf(PETSC_COMM_WORLD,"size is: %D\n",size);</i></div>
<div><i>int nstindex_f=ns*tindex_f;</i></div><div><i>ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&nstindex_f,PETSC_NULL);CHKERRQ(ierr);</i></div><div><i><br></i></div><div><i>.</i></div><div><i>.</i></div><div>
<i>.</i></div><div><i><br></i></div><div><i><br></i></div><div><i>//<span style="white-space:pre-wrap"> </span>Solution for K:</i></div><div><i>//<span style="white-space:pre-wrap"> </span>--------------</i></div>
<div><i><br></i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&nstindex_f,PETSC_NULL);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = PetscOptionsGetBool(PETSC_NULL,"-nonzero_guess",&nonzeroguess,PETSC_NULL); CHKERRQ(ierr);</i></div><div><i><br></i></div><div><i><br>
</i></div><div><i>//<span style="white-space:pre-wrap"> </span>KSPMonitorSet(ksp,MyKSPMonitor,PETSC_NULL,0);</i></div><div><i><br></i></div><div><i><br></i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = VecSetSizes(x,PETSC_DECIDE,ns*tindex_f);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = VecSetFromOptions(x);CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> </i></div><div><i><br></i></div><div><i>//<span style="white-space:pre-wrap"> </span>ierr = KSPMonitorSet(ksp,PetscInt n,PetscReal rnorm, void *);<span style="white-space:pre-wrap"> </span></i></div>
<div><i>//<span style="white-space:pre-wrap"> </span>ierr = monitor(ksp,Int it,PetscReal rnorm,);</i></div><div><i><span style="white-space:pre-wrap"> </span> </i></div><div><i><br></i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = MatTranspose(Pzz,MAT_INITIAL_MATRIX,&PzzT);</i></div><div><i><span style="white-space:pre-wrap"> </span> MatAssemblyBegin(PzzT,MAT_FINAL_ASSEMBLY);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> MatAssemblyEnd(PzzT,MAT_FINAL_ASSEMBLY);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = KSPSetOperators(ksp,PzzT,PzzT,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);</i></div>
<div><i><br></i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = PCSetType(pc,PCASM);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = KSPSetTolerances(ksp,1e-6,PETSC_DEFAULT,PETSC_DEFAULT,1e+9);CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> if (nonzeroguess)</i></div><div><i><span style="white-space:pre-wrap"> </span> {</i></div><div><i> <span style="white-space:pre-wrap"> </span>PetscScalar p = 0.5;</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = VecSet(x,p);CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> }</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = MatTranspose(Pzy,MAT_INITIAL_MATRIX,&PzyT);</i></div><div>
<i><span style="white-space:pre-wrap"> </span> MatAssemblyBegin(PzyT,MAT_FINAL_ASSEMBLY);</i></div><div><i><span style="white-space:pre-wrap"> </span> MatAssemblyEnd(PzyT,MAT_FINAL_ASSEMBLY);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = MatCreate(PETSC_COMM_WORLD,&KT);CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = MatSetSizes(KT,PETSC_DECIDE,PETSC_DECIDE,ns*tindex_f,2);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = MatSetFromOptions(KT);CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> int ctr=0;</i></div><div>
<i><span style="white-space:pre-wrap"> </span> for (ctr=0; ctr<2; ctr++)</i></div><div><i><span style="white-space:pre-wrap"> </span> {</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = VecCreate(PETSC_COMM_WORLD,&PzyTvec);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = VecSetSizes(PzyTvec,PETSC_DECIDE,ns*tindex_f);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = VecSetFromOptions(PzyTvec);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = MatGetColumnVector(PzyT,PzyTvec,ctr);</i></div><div><i><br></i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = KSPSolve(ksp,PzyTvec,x);CHKERRQ(ierr);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);</i></div><div><i> <span style="white-space:pre-wrap"> </span> ierr = KSPGetResidualNorm(ksp,& norm);<span style="white-space:pre-wrap"> </span> </i></div>
<div><i><span style="white-space:pre-wrap"> </span> printf("residual norm is=%e\n",norm);</i></div><div><i> <span style="white-space:pre-wrap"> </span> printf("iteration number is=%d\n",its);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);</i></div><div><i><span style="white-space:pre-wrap"> </span> </i></div>
<div><i>//<span style="white-space:pre-wrap"> </span>printf("%d\n", ctr);</i></div><div><i><span style="white-space:pre-wrap"> </span> for (i=0; i<ns*tindex_f; i++)</i></div>
<div><i><span style="white-space:pre-wrap"> </span> {</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = VecGetValues(x,1,&i,&xvec[i]);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = MatSetValues(KT,1,&i,1,&ctr,&xvec[i],INSERT_VALUES);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> }</i></div><div><i><span style="white-space:pre-wrap"> </span> }</i></div><div><i><span style="white-space:pre-wrap"> </span> MatAssemblyBegin(KT,MAT_FINAL_ASSEMBLY);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> MatAssemblyEnd(KT,MAT_FINAL_ASSEMBLY);</i></div><div><i><span style="white-space:pre-wrap"> </span> MatTranspose(KT,MAT_INITIAL_MATRIX,&K);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);</i></div><div><i><span style="white-space:pre-wrap"> </span> MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);</i></div>
<div><i><br></i></div><div><i><span style="white-space:pre-wrap"> </span> MatPtAP(Pzz,KT,MAT_INITIAL_MATRIX,1.0,&PK);</i></div><div><i>//MatCopy(Pzz,PzzTMP,SAME_NONZERO_PATTERN);</i></div><div><i> ierr = MatAXPY(Pmat,-1,PK,DIFFERENT_NONZERO_PATTERN);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> MatView(Pmat,PETSC_VIEWER_STDOUT_WORLD);</i></div><div><i><br></i></div><div><i> ierr = MatAXPY(Ym,-1,Mean,DIFFERENT_NONZERO_PATTERN);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> MatAssemblyBegin(Ym,MAT_FINAL_ASSEMBLY);</i></div><div><i> MatAssemblyEnd(Ym,MAT_FINAL_ASSEMBLY);</i></div><div><i><br></i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = MatMatMult(K,Ym,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&KY);</i></div>
<div><i><span style="white-space:pre-wrap"> </span> ierr = MatAssemblyBegin(KY,MAT_FINAL_ASSEMBLY);</i></div><div><i><span style="white-space:pre-wrap"> </span> ierr = MatAssemblyEnd(KY,MAT_FINAL_ASSEMBLY);</i></div>
<div><i> <span style="white-space:pre-wrap"> </span> ierr = MatAXPY(XYmat,1,KY,DIFFERENT_NONZERO_PATTERN);</i></div><div><i><span style="white-space:pre-wrap"> </span> MatView(XYmat,PETSC_VIEWER_STDOUT_WORLD);</i></div>
<div><i><br></i></div><div><i>ierr = MatDestroy(&Mean);</i></div><div><i><br></i></div><div><i>ierr=PetscFinalize();CHKERRQ(ierr);</i></div><div><i><br></i></div><div><i> return 0;</i></div><div><i><br></i></div><div>
<i> }</i></div></div><div><i><br></i></div><div><br></div><div><br></div><div>I am using the following commands to compile and run the code:</div><div><br></div><div><div><i>mpicc pcasm.c -o pcasm -l petsc</i></div></div>
<div><i>mpiexec -np 2 ./pcasm</i></div><div><br></div><div>But when I run the code I get the following error message which I don't know what it means.</div><div><br></div><div><br></div><div><div>[0]PETSC ERROR: --------------------- Error Message ------------------------------------</div>
<div>[0]PETSC ERROR: Floating point exception!</div><div>[0]PETSC ERROR: Infinite or not-a-number generated in norm!</div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div>
<div>[0]PETSC ERROR: Petsc Release Version 3.2.0, Patch 3, Fri Sep 30 10:28:33 CDT 2011</div><div>[0]PETSC ERROR: See docs/changes/index.html for recent updates.</div><div>[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.</div>
<div>[0]PETSC ERROR: See docs/index.html for manual pages.</div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: ./pcasm on a linux-imp named <a href="http://k07n14.ccr.buffalo.edu" target="_blank">k07n14.ccr.buffalo.edu</a> by rm93 Wed Dec 21 15:05:33 2011</div>
<div>[0]PETSC ERROR: Libraries linked from /util/petsc/petsc-3.2-p3/linux-impi-mkl/lib</div><div>[0]PETSC ERROR: Configure run at Fri Oct 21 08:36:23 2011</div><div>[0]PETSC ERROR: Configure options --CC=/util/intel/impi/<a href="http://4.0.3.008/intel64/bin/mpiicc" target="_blank">4.0.3.008/intel64/bin/mpiicc</a> --FC=/util/intel/impi/<a href="http://4.0.3.008/intel64/bin/mpiifort" target="_blank">4.0.3.008/intel64/bin/mpiifort</a> --CXX=/util/intel/impi/<a href="http://4.0.3.008/intel64/bin/mpiicpc" target="_blank">4.0.3.008/intel64/bin/mpiicpc</a> --with-blas-lapack-dir=/util/intel/composer_xe_2011_sp1/mkl/lib/intel64 --download-hypre=1 --with-debugging=0 -PETSC_ARCH=linux-impi-mkl --with-shared-libraries=1</div>
<div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: VecNorm() line 167 in src/vec/vec/interface/rvector.c</div><div>[0]PETSC ERROR: VecNormalize() line 261 in src/vec/vec/interface/rvector.c</div>
<div>[0]PETSC ERROR: GMREScycle() line 128 in src/ksp/ksp/impls/gmres/gmres.c</div><div>[0]PETSC ERROR: KSPSolve_GMRES() line 231 in src/ksp/ksp/impls/gmres/gmres.c</div><div>[0]PETSC ERROR: KSPSolve() line 423 in src/ksp/ksp/interface/itfunc.c</div>
<div>[0]PETSC ERROR: main() line 487 in "unknowndirectory/"pcasm.c</div><div>application called MPI_Abort(MPI_COMM_WORLD, 72) - process 0</div><div>rank 0 in job 2 k07n14.ccr.buffalo.edu_51735 caused collective abort of all ranks</div>
<div> exit status of rank 0: return code 72</div></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div>Could anyone help me on this problem? Thanks in advance for any help.</div>
<div><div><br></div><div><br></div><div><br></div></div>
</blockquote></div></div></div><span><font color="#888888"><br><br clear="all"><span><font color="#888888"><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<br>
</font></span></font></span></blockquote></div><br></div></div>
</blockquote></div></div></div><div class="HOEnZb"><div class="h5"><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<br>
</div></div></blockquote></div><br></div></div>