<div dir="ltr"><div dir="ltr">On Mon, Aug 30, 2021 at 4:42 PM Stefano Zampini <<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">What is the error you are getting from the KSP? Default solver in parallel in BlockJacobi+ILU which does not work for non-square matrices. You do not need to call PCSetFromOptions on the pc. Just call KSPSetFromOptions and run with -pc_type none<br></blockquote><div><br></div><div>I have a more basic question. The idea of GMRES is the following:</div><div><br></div><div> We build a space, {r, A r, A^2 r, ...} for the solution</div><div><br></div><div>This means that r and Ar must have compatible layouts. It does not sound like this is the case for you.</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:1px solid rgb(204,204,204);padding-left:1ex">
> On Aug 30, 2021, at 9:17 PM, Olivier Jamond <<a href="mailto:olivier.jamond@cea.fr" target="_blank">olivier.jamond@cea.fr</a>> wrote:<br>
> <br>
> Hello,<br>
> <br>
> I am sorry because I surely miss something, but I cannot manage to solve a problem with a MPIAIJ matrix which has non-square 'diagonal parts'.<br>
> <br>
> I copy/paste at the bottom of this message a very simple piece of code which causes me troubles. I this code I try to do 'x=KSP(A)*b' (with gmres/jacobi), but this fails whereas a matrix multiplication 'b=A*x' seems to work. Is ksp with such a matrix supposed to work (I can't find anything in the documentation about that, so I guess that it is...) ?<br>
> <br>
> Many thanks,<br>
> Olivier<br>
> <br>
> NB: this code should be launched with exactly 3 procs<br>
> <br>
> #include "petscsys.h" /* framework routines */<br>
> #include "petscvec.h" /* vectors */<br>
> #include "petscmat.h" /* matrices */<br>
> #include "petscksp.h"<br>
> <br>
> #include <vector><br>
> #include <string><br>
> #include <iostream><br>
> #include <numeric><br>
> <br>
> static char help[] = "Trying to solve a linear system on a sub-block of a matrix\n\n";<br>
> int main(int argc, char **argv)<br>
> {<br>
> MPI_Init(NULL, NULL);<br>
> PetscErrorCode ierr;<br>
> ierr = PetscInitialize(&argc, &argv, NULL, help);<br>
> CHKERRQ(ierr);<br>
> <br>
> // clang-format off<br>
> std::vector<std::vector<double>> AA = {<br>
> { 1, 2, 0, /**/ 0, 3, /**/ 0, 0, 4},<br>
> { 0, 5, 6, /**/ 7, 0, /**/ 0, 8, 0},<br>
> { 9, 0, 10, /**/11, 0, /**/ 0, 12, 0},<br>
> //---------------------------------------<br>
> {13, 0, 14, /**/15, 16, /**/17, 0, 0},<br>
> { 0, 18, 0, /**/19, 20, /**/21, 0, 0},<br>
> { 0, 0, 0, /**/22, 23, /**/ 1, 24, 0},<br>
> //--------------------------------------<br>
> {25, 26, 27, /**/ 0, 0, /**/28, 29, 0},<br>
> {30, 0, 0, /**/31, 32, /**/33, 0, 34},<br>
> };<br>
> <br>
> <br>
> std::vector<double> bb = {1.,<br>
> 1.,<br>
> 1.,<br>
> //<br>
> 1.,<br>
> 1.,<br>
> 1.,<br>
> //<br>
> 1.,<br>
> 1.};<br>
> <br>
> <br>
> std::vector<int> nDofsRow = {3, 3, 2};<br>
> std::vector<int> nDofsCol = {3, 2, 3};<br>
> // clang-format on<br>
> <br>
> int NDofs = std::accumulate(nDofsRow.begin(), nDofsRow.end(), 0);<br>
> <br>
> int pRank, nProc;<br>
> MPI_Comm_rank(PETSC_COMM_WORLD, &pRank);<br>
> MPI_Comm_size(PETSC_COMM_WORLD, &nProc);<br>
> <br>
> if (nProc != 3) {<br>
> std::cerr << "THIS TEST MUST BE LAUNCHED WITH EXACTLY 3 PROCS\n";<br>
> abort();<br>
> }<br>
> <br>
> Mat A;<br>
> MatCreate(PETSC_COMM_WORLD, &A);<br>
> MatSetType(A, MATMPIAIJ);<br>
> MatSetSizes(A, nDofsRow[pRank], nDofsCol[pRank], PETSC_DETERMINE, PETSC_DETERMINE);<br>
> MatMPIAIJSetPreallocation(A, NDofs, NULL, NDofs, NULL);<br>
> <br>
> Vec b;<br>
> VecCreate(PETSC_COMM_WORLD, &b);<br>
> VecSetType(b, VECMPI);<br>
> VecSetSizes(b, nDofsRow[pRank], PETSC_DECIDE);<br>
> <br>
> if (pRank == 0) {<br>
> for (int i = 0; i < NDofs; ++i) {<br>
> for (int j = 0; j < NDofs; ++j) {<br>
> if (AA[i][j] != 0.) {<br>
> MatSetValue(A, i, j, AA[i][j], ADD_VALUES);<br>
> }<br>
> }<br>
> VecSetValue(b, i, bb[i], ADD_VALUES);<br>
> }<br>
> }<br>
> <br>
> MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);<br>
> MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);<br>
> VecAssemblyBegin(b);<br>
> VecAssemblyEnd(b);<br>
> <br>
> PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_DENSE);<br>
> MatView(A, PETSC_VIEWER_STDOUT_WORLD);<br>
> VecView(b, PETSC_VIEWER_STDOUT_WORLD);<br>
> <br>
> KSP ksp;<br>
> KSPCreate(PETSC_COMM_WORLD, &ksp);<br>
> KSPSetOperators(ksp, A, A);<br>
> KSPSetFromOptions(ksp);<br>
> <br>
> PC pc;<br>
> KSPGetPC(ksp, &pc);<br>
> PCSetFromOptions(pc);<br>
> <br>
> Vec x;<br>
> MatCreateVecs(A, &x, NULL);<br>
> ierr = KSPSolve(ksp, b, x); // this fails<br>
> MatMult(A, x, b); // whereas the seems to be ok...<br>
> <br>
> VecView(x, PETSC_VIEWER_STDOUT_WORLD);<br>
> <br>
> MPI_Finalize();<br>
> <br>
> return 0;<br>
> }<br>
> <br>
<br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>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><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>