Hi Jed<div><br></div><div>After probing why the following code snippit isn't run:</div><div><div style="color:rgb(80,0,80);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)"><br class="Apple-interchange-newline">
if (mg->galerkin == 1) {</div><div style="color:rgb(80,0,80);font-family:arial,sans-serif;font-size:13px;background-color:rgb(255,255,255)"><div> Mat B;</div><div> /* currently only handle case where mat and pmat are the same on coarser levels */</div>
<div> ierr = KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB,&uflag);CHKERRQ(ierr);</div><div> if (!pc->setupcalled) {</div><div> for (i=n-2; i>-1; i--) {</div><div> ierr = MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);CHKERRQ(ierr);</div>
<div> ierr = KSPSetOperators(mglevels[i]->smoothd,B,B,uflag);CHKERRQ(ierr);</div><div> if (i != n-2) {ierr = PetscObjectDereference((PetscObject)dB);CHKERRQ(ierr);}</div><div> dB = B;</div><div> }</div>
<div> if (n > 1) {ierr = PetscObjectDereference((PetscObject)dB);CHKERRQ(ierr);}</div><div> } else {</div><div> for (i=n-2; i>-1; i--) {</div><div> ierr = KSPGetOperators(mglevels[i]->smoothd,PETSC_NULL,&B,PETSC_NULL);CHKERRQ(ierr);</div>
<div> ierr = MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);CHKERRQ(ierr);</div><div> ierr = KSPSetOperators(mglevels[i]->smoothd,B,B,uflag);CHKERRQ(ierr);</div><div> dB = B;</div>
<div> }</div><div> }</div><div><br></div></div>It is because after I use PCMGSetGalerkin(), mg->galerkin has a value of -1 not 1. From what I can tell the value of PETSC_TRUE is architecture/implementation dependant. I can't see where/how the value of mg->galerkin obtains a value of 1 and runs the code above. </div>
<div><br></div><div>As for the other issue, after converting my baij matrix to aij, running MatRART() is ok. </div><div><br></div><div>So we just have to figure out why calling</div><div><br></div><div><br></div><div><div>
call PCMGSetGalerkin(pc, PETSC_TRUE, PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div></div><div><br></div><div>sets mg->galerkin to a value of -1 and the code above is expecting a value of 1.</div>
<div><br></div><div>Gaetan</div><div><br><br><br><br><div class="gmail_quote">On Sat, Sep 1, 2012 at 5:34 PM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><div>On Sat, Sep 1, 2012 at 4:20 PM, Gaetan Kenway <span dir="ltr"><<a href="mailto:kenway@utias.utoronto.ca" target="_blank">kenway@utias.utoronto.ca</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>Sorry about the last email...I hit sent it too early. </div><div><div><br></div><div>The full calling sequence is:</div><div><br></div></div><div><div><div><div><br> call KSPCreate(SUMB_COMM_WORLD, ksp, PETScIerr)</div>
<div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><br></div><div> useAD = .False.</div><div> useTranspose = .True.</div><div> usePC = .True.</div><div> call setupStateResidualMatrix(drdwpret,useAD,usePC,useTranspose)</div>
<div><br></div></div><div> ! dRdwT has already been created</div><div><div><div> call KSPSetOperators(ksp, dRdWT, dRdwpreT, &</div><div> DIFFERENT_NONZERO_PATTERN, PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div>
<div><br></div><div> call KSPSetFromOptions(ksp, PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><br></div><div> call KSPSetType(ksp, KSPFGMRES, PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div>
<div><br></div><div> call KSPGMRESSetRestart(ksp, adjRestart, PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><br></div><div> call KSPGMRESSetCGSRefinementType(ksp, &</div><div> KSP_GMRES_CGS_REFINE_IFNEEDED, PETScIerr)</div>
<div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><br></div><div> call kspgetpc(ksp, pc, PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><br></div><div> call pcsettype(pc, PCMG, PETScIerr)</div>
<div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><br></div><div> call PetscOptionsSetValue('-pc_mg_levels','2',PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div></div></div></div>
<div><div><div><div><br>
</div><div> ! Set command line options </div>
<div> call PCSetFromOptions(pc, PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><div><br></div><div> call PCMGSetGalerkin(pc, PETSC_TRUE, PETScIerr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div>
<div><br></div></div><div> ! Create the restriction operator between finest and one below: </div>
<div> call createRestrictOperator(RL1, 1)</div><div><br></div><div> call PCMGSetRestriction(pc, 1, RL1, PETScierr)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><br></div><div> call PCSetup(pc, PETScierr)</div>
<div> call EChk(PETScIerr,__FILE__,__LINE__)</div></div><div><div></div></div></div></div></div><div><br></div><div> The error occurs in PCSetup() which calls PCSetup_MG() and the issue is in KSPGetVecs() on line 644 in mg.c. Just to clarify, this isn't a memory problem/segfault.</div>
</blockquote><div><br></div></div></div><div>The output below is a SEGV.</div><div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div> PETSc returns control cleanly and I my code exits when the error is detected. Therefore a backtrace or valgrind doesn't yield any additional information. ( I Ran both).</div>
</blockquote><div><br></div></div><div>You can still set a breakpoint in PetscError (-start_in_debugger / -on_error_attach_debugger does this automatically). Then you can examine the stack. </div><div><div> </div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div> The terminal output is:</div>
<div><br></div><div><div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range</div>
<div>[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger</div><div>[0]PETSC ERROR: or see <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC</a> ERROR: or try <a href="http://valgrind.org" target="_blank">http://valgrind.org</a> on GNU/linux and Apple Mac OS X to find memory corruption errors</div>
<div>[0]PETSC ERROR: likely location of problem given in stack below</div><div>[0]PETSC ERROR: --------------------- Stack Frames ------------------------------------</div></div><div><div> =================================================================</div>
<div>PETSc or MPI Error. Error Code 73. Detected on Proc 0</div></div><div>Error at line: 130 in file: setupPETScKsp.F90</div><div> =================================================================</div><div>
<div>[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,</div>
<div>[0]PETSC ERROR: INSTEAD the line number of the start of the function</div><div>[0]PETSC ERROR: is given.</div><div>[0]PETSC ERROR: [0] MatGetVecs line 8142 /nfs/mica/home/kenway/Downloads/petsc-3.3/src/mat/interface/matrix.c</div>
<div>[0]PETSC ERROR: [0] KSPGetVecs line 774 /nfs/mica/home/kenway/Downloads/petsc-3.3/src/ksp/ksp/interface/iterativ.c</div><div>[0]PETSC ERROR: [0] PCSetUp_MG line 508 /nfs/mica/home/kenway/Downloads/petsc-3.3/src/ksp/pc/impls/mg/mg.c</div>
<div>[0]PETSC ERROR: [0] PCSetUp line 810 /nfs/mica/home/kenway/Downloads/petsc-3.3/src/ksp/pc/interface/precon.c</div></div></div><div><br></div>I had a look at ex42.c but it uses a DM which I am not doing.</blockquote>
<div><br></div></div><div>It uses the DM for convenience to build the operators, but it does *not* tell the KSP about the DM. That is, it does not call KSPSetDM() or PCSetDM(). So the code in mg.c does not know about the DM.</div>
<div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> From the example there is the following:<div>
<br></div><div><div>ierr = PCMGSetLevels(pc,nlevels,PETSC_NULL);CHKERRQ(ierr);</div><div> ierr = PCMGSetType(pc,PC_MG_MULTIPLICATIVE);CHKERRQ(ierr);</div><div> ierr = PCMGSetGalerkin(pc,PETSC_TRUE);CHKERRQ(ierr);</div>
<div>
<br></div><div> for (k=1; k<nlevels; k++){</div><div> ierr = DMCreateInterpolation(da_list[k-1],da_list[k],&R,PETSC_NULL);CHKERRQ(ierr);</div><div> ierr = PCMGSetInterpolation(pc,k,R);CHKERRQ(ierr);</div><div>
ierr = MatDestroy(&R);CHKERRQ(ierr);</div><div> }</div><div><br></div><div>Is there anything special about "R" in this case?</div></div></blockquote><div><br></div></div><div>Nothing special, it is just an MAIJ matrix that performs the interpolation.</div>
<div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div> As I mentioned before, I know what the issue is, the coarse grid operators are not set at all. From what I can tell, they should be set in the routine PCSetup_MG in mg.c, but the code is never run unless pc->dm is True. </div>
</div></blockquote><div><br></div></div><div>What code is not run unless pc->dm?</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 guess I'll just create the operators myself. Even if I modify the code to make it enter the if statement at line 577</div>
</div></blockquote><div><br></div></div><div>You called PCMGSetGalerkin() so the code below should run, constructing the coarse levels. You could try setting a breakpoint in there to see why it's not running.</div><div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>
<div><br></div><div><div> if (mg->galerkin == 1) {</div></div><div><div> Mat B;</div><div> /* currently only handle case where mat and pmat are the same on coarser levels */</div><div> ierr = KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB,&uflag);CHKERRQ(ierr);</div>
<div> if (!pc->setupcalled) {</div><div> for (i=n-2; i>-1; i--) {</div><div> ierr = MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);CHKERRQ(ierr);</div><div> ierr = KSPSetOperators(mglevels[i]->smoothd,B,B,uflag);CHKERRQ(ierr);</div>
<div> if (i != n-2) {ierr = PetscObjectDereference((PetscObject)dB);CHKERRQ(ierr);}</div><div> dB = B;</div><div> }</div><div> if (n > 1) {ierr = PetscObjectDereference((PetscObject)dB);CHKERRQ(ierr);}</div>
<div> } else {</div><div> for (i=n-2; i>-1; i--) {</div><div> ierr = KSPGetOperators(mglevels[i]->smoothd,PETSC_NULL,&B,PETSC_NULL);CHKERRQ(ierr);</div><div> ierr = MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);CHKERRQ(ierr);</div>
<div> ierr = KSPSetOperators(mglevels[i]->smoothd,B,B,uflag);CHKERRQ(ierr);</div><div> dB = B;</div><div> }</div><div> }</div></div><div><br></div><div>which seems like where it should be done....the MatPTAP will not run.</div>
</div></blockquote><div><br></div></div><div>What do you mean "will not run"?</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'm using a mpibaij matrix for 'dRdwPreT' and my restriction matrix 'RL1' is a mpiaij matrix. When I try to run MatRART myself with the code:</div>
<div><br></div><div><div> call MatRART(dRdwT, RL1, MAT_INITIAL_MATRIX, 1.0, newMat, PETScIERR)</div><div> call EChk(PETScIerr,__FILE__,__LINE__)</div><div><br></div><div>I receive error code 56: /* no support for requested operation */ (this is a clean exit from PETSc)</div>
</div></div></blockquote><div><br></div></div><div>MatRARt (or MatPtAP, for that matter) is not implemented for BAIJ with an AIJ projection. The AIJ projection will lose the block structure anyway, so just use AIJ for both.</div>
<div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div>
</div><div><br></div><div>In conclusion, I have two issues: 1) The coarse grid operators are not being set automatically in mg.c. 2) Using MatRART() with a mpibaij matrix and a mpiaij matrix is unsupported. </div><div><br>
</div><div>Is there a PCMG example that doesn't use a DM to generate the restriction/prolongation/matrix information I can have a look at?</div></div></blockquote></div></div><br>
</blockquote></div><br>
</div>