I'm attempting to integrate PETSc into an existing simulation and I have successfully gotten PETSc to handle the first matrix inversion, but I'm having trouble setting it up properly for subsequent inversions. It currently does the inversions without complaining, but gives incorrect output. I believe I need to 'reset' the matrix being inverted. I am inverting the same matrix each time, but the values and nonzero structure change. I create and set up the matrix, ksp, and pc once, then call a routine that fills the matrix and calls KSPSetOperators() before using MatMatSolve again. The routine is called multiple times.<div>
<br></div><div>I am including the actual code, but here is the outline of what I'm currently trying:<div><br></div><div><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo;color:#138607">//Setup hamiltonian</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = MatCreate(PETSC_COMM_WORLD, &hamiltonian);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = MatSetType(hamiltonian, MATAIJ);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = MatSetSizes(hamiltonian, PETSC_DECIDE, PETSC_DECIDE, rows, columns);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = MatAssemblyBegin(hamiltonian, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = MatAssemblyEnd(hamiltonian, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo"><br></p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">
</p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo;color:#138607">//Setup PETSc Solver</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = KSPSetOperators(ksp, hamiltonian, hamiltonian, SAME_NONZERO_PATTERN);CHKERRQ(ierr); </p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);</p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo"></p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">
ierr = test_invert(hamiltonian, inverseHamiltonian, identityMat, pc, ksp, ...);CHKERRQ(ierr);</p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo"><br></p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">
<br></p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo"><font face="arial, helvetica, sans-serif"><span style="font-size:small">test_invert (looped</span></font><span style="font-family:arial, helvetica, sans-serif;font-size:small"> part):</span></p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo"><font face="arial, helvetica, sans-serif"><span style="font-size:small"><br></span></font></p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">
<font face="arial, helvetica, sans-serif"><span style="font-size:small"><span style="font-family:Menlo;font-size:11px;color:rgb(19, 134, 7)">//Several MatSetValue() calls</span></span></font></p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo"></p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo;color:#138607">//Setup PETSc Hamiltonian</p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">
ierr = MatAssemblyBegin(hamiltonian, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = MatAssemblyEnd(hamiltonian, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo;min-height:13.0px"><span style="white-space:pre-wrap"> </span></p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo;color:#138607">//Setup PETSc Solver</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = KSPSetOperators(ksp, hamiltonian, hamiltonian, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = KSPSetUp(ksp);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo;min-height:13.0px"><span style="white-space:pre-wrap"> </span></p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo;color:#138607">//Factor and solve</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = PCFactorGetMatrix(pc, &hamiltonian);CHKERRQ(ierr);</p>
<p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">ierr = MatMatSolve(hamiltonian, identityMat, inverseHamiltonian);CHKERRQ(ierr);</p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo"><br></p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo">
<br></p><p style="margin:0.0px 0.0px 0.0px 0.0px;font:11.0px Menlo"><font face="arial, helvetica, sans-serif"><span style="font-size:small">What needs to be done to the matrix to clear it out and use it again? MatZeroEntries doesn't want to work on an unfactored matrix, plus the nonzero structure changes. Do I need to reset the inverseHamiltonian as well?</span></font></p>
<p></p><p></p></div></div>