<div dir="ltr">Hey, <div>i tested the fix today and ran the same code as above. It is now working fine and my results are in agreement with my analytical solutions. Although i have to tranpose the computed entries afterwards. This should be the correct behaviour or ? <br> ierr = MatTranspose(spRHST, MAT_INPLACE_MATRIX, &spRHST);CHKERRQ(ierr);</div><div><br><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Am Mo., 16. Dez. 2019 um 22:57 Uhr schrieb Zhang, Hong <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>>:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
<div dir="ltr">
<div dir="ltr">Jan:<br>
</div>
<div>You found a bug in our interface. I pushed a fix <a href="https://gitlab.com/petsc/petsc/commit/6fa4709e1488f7c58b0a4ecbdf34b75ae4ac64d2" target="_blank">https://gitlab.com/petsc/petsc/commit/6fa4709e1488f7c58b0a4ecbdf34b75ae4ac64d2</a></div>
<div>in branch hzhang/fix-mumps-GetInverse. </div>
<div>You may give it a try by 'git checkout hzhang/fix-mumps-GetInverse' at PETSC_DIR and rebuild petsc.</div>
<div>Let me know if it fixes the issue.</div>
<div>Hong</div>
<div><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">
<div dir="ltr">Hello, everybody, <br>
I am using PETSc version 3.12 together with MUMPS to calculate a part of an inverse from a matrix A. For this I used the example mat/examples/tests/ex214.c as suggested here in the forum. For test purposes I read a matrix with the dimension 10x10, which only
has entries on the main diagonal, because I know the inverse of it. <br>
But now I have a problem using the MatMumpsGetInverse() function. I get the error message:<br>
Error reported by MUMPS in solve phase: INFOG(1)=-47 INFO(2)=1 ACcording to the MUMPS Docu this message tells me that I would ignored the constraint NRHS= N. I checked the shape of spRHST but they appear to be correct. Does anyone see the error?
<br>
My code is appended below.
<div>static char help[] ="Compute a part of the inverse of a sparse matrix. This code requires that PETSc was configured with MUMPS since we are dealing with large matrices \<br>
and therefore use a parallel LU factorization. We compute the inverse by solving the equation A*X=RHS. Where A is our Matrix, X is the inverse and RHS is the identity matrix.\<br>
Note that the number of columns nrhs of X can be chosen smaller than the number of columns N in A. Therefore only a part of the inverse is computed in X. \n \<br>
In this code we use a sparse representation of the RHS matrix in MUMPS in csr format. Computation of selected entries in inv(A) is done using MatMumpsGetInverse. \n \<br>
Input parameters: \n\<br>
-fin <input_file> : file to load \n \<br>
-fout <input_file> : file to load \n \<br>
-nrhs <numberofcolumns> : Number of columns to compute \n \<br>
-displ <Bool>: Print matrices to terminal \n\<br>
Example usage: \n \<br>
mpiexec -np 2 ./compute_inverse_sparse_rhs -fin ../../convert_to_binary_petsc_matrix/identity_matrix_prefactor3_ncols10 -nrhs 5 -displ";<br>
<br>
#include <stdio.h><br>
#include <petscmat.h><br>
#include <petscviewer.h> <br>
<br>
int main(int argc, char **args){<br>
PetscErrorCode ierr; // Datatype used for return error code<br>
PetscMPIInt size,rank; // Datatype used to represent 'int' parameters to MPI functions.<br>
#if defined(PETSC_HAVE_MUMPS)<br>
Mat A,F,spRHST; // Abstract PETSc matrix object used to manage all linear operators in PETSc<br>
PetscViewer fd; // Abstract PETSc object that helps view (in ASCII, binary, graphically etc) other PETSc objects<br>
PetscBool flg1,flg2; // Logical variable. Actually an int in C.<br>
PetscBool displ=PETSC_FALSE; // Display matrices <br>
PetscInt M,N,m,n,rstart,rend,nrhs,i; // PETSc type that represents an integer, used primarily to represent size of arrays and indexing into arrays.<br>
PetscScalar v; // PETSc type that represents either a double precision real number,...<br>
char inputfile[1][PETSC_MAX_PATH_LEN]; // Input file name <br>
// char outputfile[1][PETSC_MAX_PATH_LEN]; // Outputfile file name <br>
#endif<br>
<br>
// Initializes PETSc and MPI. Get size and rank of MPI.<br>
ierr = PetscInitialize(&argc, &args, (char*)0, help);if (ierr){return ierr;}<br>
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);<br>
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);<br>
<br>
//Check if PETSc was configured with MUMPS. If not print error message and exit
<br>
#if !defined(PETSC_HAVE_MUMPS)<br>
if (!=rank){ierr = PetscPrintf(PETSC_COMM_SELF, "This code requires MUMPS, exit...\n");CHKERRQ(ierr);<br>
ierr = PetscFinalize();<br>
return ierr;<br>
}<br>
#else<br>
<br>
// Check if displ is set. If True the matrices are printed to the terminal<br>
ierr = PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL);CHKERRQ(ierr);<br>
<br>
// Load matrix A from file <br>
ierr = PetscOptionsGetString(NULL, NULL, "-fin" ,inputfile[0], PETSC_MAX_PATH_LEN, &flg1);CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "Load matrix in: %s \n", inputfile[0]);CHKERRQ(ierr);<br>
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, inputfile[0], FILE_MODE_READ, &fd);CHKERRQ(ierr);<br>
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);<br>
ierr = MatSetType(A, MATAIJ);CHKERRQ(ierr); <br>
ierr = MatLoad(A, fd);CHKERRQ(ierr);<br>
// Print matrix A <br>
if (displ){<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "\n---------------\n");CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD,"Matrix A from file:\n", nrhs);<br>
ierr = MatView(A, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "\n");CHKERRQ(ierr);<br>
}<br>
// Check if matrix is quadratic<br>
ierr = MatGetSize(A, &M, &N);CHKERRQ(ierr);<br>
ierr = MatGetLocalSize(A, &m, &n);CHKERRQ(ierr);<br>
if (M != N){<br>
//Macro that is called when an error has been detected.<br>
SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Detected a rectangular matrix: (%d, %d)", M, N);<br>
}<br>
ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "---------------\n");CHKERRQ(ierr);<br>
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Ownership ranges for Matrix A, rank: %i, size: %i, rstart: %i, rend: %i \n", rank, size, rstart, rend);CHKERRQ(ierr);<br>
ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "---------------\n");CHKERRQ(ierr);<br>
<br>
// Set the number of columns of the inverse to be computed. <br>
nrhs = N;<br>
ierr = PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, &flg2);CHKERRQ(ierr);<br>
<br>
// Create SpRHST for inv(A) with sparse RHS stored in the host.<br>
// PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix,<br>
// thus user must create spRHST=spRHS^T and call MatMatTransposeSolve()<br>
// User must create B^T in sparse compressed row format on the host processor and call MatMatTransposeSolve() to implement MUMPS' MatMatSolve().<br>
// MUMPS requires nrhs = N <br>
ierr = MatCreate(PETSC_COMM_WORLD, &spRHST);CHKERRQ(ierr);<br>
if (!rank){<br>
ierr = MatSetSizes(spRHST,N,M,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);<br>
}<br>
else{<br>
ierr = MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);<br>
}<br>
ierr = MatSetType(spRHST,MATAIJ);CHKERRQ(ierr);<br>
ierr = MatSetFromOptions(spRHST);CHKERRQ(ierr);<br>
ierr = MatSetUp(spRHST);CHKERRQ(ierr);<br>
if (!rank){<br>
v = 1.0;<br>
for(i=0;i<nrhs;i++){<br>
ierr = MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);<br>
}<br>
}<br>
ierr = MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(spRHST, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br>
<br>
// Print matrix spRHST <br>
if (displ){<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "\n---------------\n");CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD,"Matrix spRHST:\n", nrhs);<br>
ierr = MatView(spRHST, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "\n");CHKERRQ(ierr);<br>
}<br>
<br>
// Print information<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "\nCompute %i columns of the inverse using LU-factorization in MUMPS!\n", nrhs);<br>
<br>
// Factorize the Matrix using a parallel LU factorization in MUMPS<br>
ierr = MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F);CHKERRQ(ierr);<br>
ierr = MatLUFactorSymbolic(F, A, NULL, NULL, NULL);CHKERRQ(ierr);<br>
ierr = MatLUFactorNumeric(F, A, NULL);CHKERRQ(ierr);<br>
<br>
// Create spRHS <br>
Mat spRHS = NULL;<br>
<br>
// Create spRHS = spRHS^T. Two matrices that share internal matrix data structure.
<br>
// Creates a new matrix object that behaves like A'.<br>
ierr = MatCreateTranspose(spRHST,&spRHS);CHKERRQ(ierr);<br>
<br>
// Get user-specified set of entries in inverse of A<br>
ierr = MatMumpsGetInverse(F,spRHS);CHKERRQ(ierr);<br>
<br>
if (displ){<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "\n---------------\n");CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD,"First %D columns of inv(A) with sparse RHS:\n", nrhs);<br>
ierr = MatView(spRHS,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br>
ierr = PetscPrintf(PETSC_COMM_WORLD, "---------------\n");CHKERRQ(ierr);<br>
}<br>
<br>
// Free data structures<br>
ierr = MatDestroy(&A);CHKERRQ(ierr);<br>
ierr = MatDestroy(&spRHS);CHKERRQ(ierr);<br>
ierr = MatDestroy(&spRHST);CHKERRQ(ierr);<br>
ierr = PetscFinalize();<br>
return ierr;<br>
<br>
#endif<br>
}<br>
<br>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</blockquote></div>