<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
</head>
<body>
<div dir="ltr">
<div dir="ltr">Hello,
<div><br>
</div>
<div>For some small matrices, I'd like to use Petsc to perform direct LU factorization on a sequential dense or AIJ matrix, and then use the factored matrix later on via MatSolve. This occurs multiple times in my code, and each factored matrix is in turn used
multiple times.</div>
<div><br>
</div>
<div>I tried to wrap this factorization process in a function, which should return the factored matrix, as follows:</div>
<div><br>
</div>
<div>
<div>void MatFactorize_Petsc(Mat &mat, Mat &mat_factored)</div>
<div>{</div>
<div><span style="white-space:pre-wrap"></span>PC pc;</div>
<div><span style="white-space:pre-wrap"></span>PCCreate(MPI_COMM_SELF, &pc);</div>
<div><span style="white-space:pre-wrap"></span>PCSetOperators(pc, mat, mat);</div>
<div><span style="white-space:pre-wrap"></span>PCSetType(pc, PCLU);</div>
<div><br>
</div>
<div><span style="white-space:pre-wrap"></span>PCFactorSetMatSolverType(pc, MATSOLVERPETSC); // Or SuperLU</div>
<div><span style="white-space:pre-wrap"></span>PCFactorSetUpMatSolverType(pc);</div>
<div><br>
</div>
<div><span style="white-space:pre-wrap"></span>PCFactorGetMatrix(pc, &mat_factored);<span style="white-space:pre-wrap">
</span></div>
<div><span style="white-space:pre-wrap"></span>PCSetUp(pc);</div>
<div><span style="white-space:pre-wrap"></span></div>
<div><span style="white-space:pre-wrap"></span>// PCDestroy(&pc);</div>
<div><span style="white-space:pre-wrap"></span></div>
<div><span style="white-space:pre-wrap"></span>return;</div>
<div>}</div>
</div>
<div><br>
</div>
<div>The command PCDestroy causes a segmentation fault, which I think happens because retrieving the factored matrix does not increase the reference count. Looking at the Petsc source, it basically returns a pointer to pc->data. So if I want to use mat_factored
outside the function, I cannot destroy the PC object, which leads to memory leaks (as per Valgrind) even if I later call MatDestroy(&mat_factored).</div>
<div><br>
</div>
<div>I tried using a temp matrix to get the factors, and then doing MatDuplicate(temp, MAT_COPY_VALUES, &mat_factored), but MatDuplicate is not allowed on factored matrices; same for MatConvert.</div>
<div><br>
</div>
<div>Is there a way I can achieve the desired behaviour, where mat_factored is not "linked" to the pc or ksp object, keeping in mind that I'd like to be able to choose between SuperLU, Petsc and SuperLU_dist solvers?</div>
<div><br>
</div>
<div>Thanks,</div>
<div>Shash</div>
<div>
<div>
<div dir="ltr">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div dir="ltr"><br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</body>
</html>