Jan: Hey, 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. Good. I'll push the fix to the petsc release after regression tests.
Although i have to tranpose the computed entries afterwards. This should be the correct behaviour or ? ierr = MatTranspose(spRHST, MAT_INPLACE_MATRIX, &spRHST);CHKERRQ(ierr); Yes, you have to do so. The issue is, mumps uses centralized compressed COLUMN format for RHS, while pestc only supports compressed ROW format. A hack is to have user create spRHST, then call MatTranspose(spRHST, MAT_INPLACE_MATRIX, &spRHS); Note: there is no copying in this routine. spRHST and spRHS share same data structure using column-wise and row-wise orderings. If you have a better suggestion, let us know. Thanks for reporting the bug. Hong Am Mo., 16. Dez. 2019 um 22:57 Uhr schrieb Zhang, Hong <[email protected]<mailto:[email protected]>>: Jan: You found a bug in our interface. I pushed a fix https://gitlab.com/petsc/petsc/commit/6fa4709e1488f7c58b0a4ecbdf34b75ae4ac64d2 in branch hzhang/fix-mumps-GetInverse. You may give it a try by 'git checkout hzhang/fix-mumps-GetInverse' at PETSC_DIR and rebuild petsc. Let me know if it fixes the issue. Hong Hello, everybody, 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. But now I have a problem using the MatMumpsGetInverse() function. I get the error message: 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? My code is appended below. 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 \ 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.\ 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 \ 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 \ Input parameters: \n\ -fin <input_file> : file to load \n \ -fout <input_file> : file to load \n \ -nrhs <numberofcolumns> : Number of columns to compute \n \ -displ <Bool>: Print matrices to terminal \n\ Example usage: \n \ mpiexec -np 2 ./compute_inverse_sparse_rhs -fin ../../convert_to_binary_petsc_matrix/identity_matrix_prefactor3_ncols10 -nrhs 5 -displ"; #include <stdio.h> #include <petscmat.h> #include <petscviewer.h> int main(int argc, char **args){ PetscErrorCode ierr; // Datatype used for return error code PetscMPIInt size,rank; // Datatype used to represent 'int' parameters to MPI functions. #if defined(PETSC_HAVE_MUMPS) Mat A,F,spRHST; // Abstract PETSc matrix object used to manage all linear operators in PETSc PetscViewer fd; // Abstract PETSc object that helps view (in ASCII, binary, graphically etc) other PETSc objects PetscBool flg1,flg2; // Logical variable. Actually an int in C. PetscBool displ=PETSC_FALSE; // Display matrices 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. PetscScalar v; // PETSc type that represents either a double precision real number,... char inputfile[1][PETSC_MAX_PATH_LEN]; // Input file name // char outputfile[1][PETSC_MAX_PATH_LEN]; // Outputfile file name #endif // Initializes PETSc and MPI. Get size and rank of MPI. ierr = PetscInitialize(&argc, &args, (char*)0, help);if (ierr){return ierr;} ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr); ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); //Check if PETSc was configured with MUMPS. If not print error message and exit #if !defined(PETSC_HAVE_MUMPS) if (!=rank){ierr = PetscPrintf(PETSC_COMM_SELF, "This code requires MUMPS, exit...\n");CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; } #else // Check if displ is set. If True the matrices are printed to the terminal ierr = PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL);CHKERRQ(ierr); // Load matrix A from file ierr = PetscOptionsGetString(NULL, NULL, "-fin" ,inputfile[0], PETSC_MAX_PATH_LEN, &flg1);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, "Load matrix in: %s \n", inputfile[0]);CHKERRQ(ierr); ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, inputfile[0], FILE_MODE_READ, &fd);CHKERRQ(ierr); ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); ierr = MatSetType(A, MATAIJ);CHKERRQ(ierr); ierr = MatLoad(A, fd);CHKERRQ(ierr); // Print matrix A if (displ){ ierr = PetscPrintf(PETSC_COMM_WORLD, "\n---------------\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Matrix A from file:\n", nrhs); ierr = MatView(A, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, "\n");CHKERRQ(ierr); } // Check if matrix is quadratic ierr = MatGetSize(A, &M, &N);CHKERRQ(ierr); ierr = MatGetLocalSize(A, &m, &n);CHKERRQ(ierr); if (M != N){ //Macro that is called when an error has been detected. SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Detected a rectangular matrix: (%d, %d)", M, N); } ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, "---------------\n");CHKERRQ(ierr); 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); ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, "---------------\n");CHKERRQ(ierr); // Set the number of columns of the inverse to be computed. nrhs = N; ierr = PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, &flg2);CHKERRQ(ierr); // Create SpRHST for inv(A) with sparse RHS stored in the host. // PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix, // thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() // User must create B^T in sparse compressed row format on the host processor and call MatMatTransposeSolve() to implement MUMPS' MatMatSolve(). // MUMPS requires nrhs = N ierr = MatCreate(PETSC_COMM_WORLD, &spRHST);CHKERRQ(ierr); if (!rank){ ierr = MatSetSizes(spRHST,N,M,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); } else{ ierr = MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); } ierr = MatSetType(spRHST,MATAIJ);CHKERRQ(ierr); ierr = MatSetFromOptions(spRHST);CHKERRQ(ierr); ierr = MatSetUp(spRHST);CHKERRQ(ierr); if (!rank){ v = 1.0; for(i=0;i<nrhs;i++){ ierr = MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); } } ierr = MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(spRHST, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); // Print matrix spRHST if (displ){ ierr = PetscPrintf(PETSC_COMM_WORLD, "\n---------------\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"Matrix spRHST:\n", nrhs); ierr = MatView(spRHST, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, "\n");CHKERRQ(ierr); } // Print information ierr = PetscPrintf(PETSC_COMM_WORLD, "\nCompute %i columns of the inverse using LU-factorization in MUMPS!\n", nrhs); // Factorize the Matrix using a parallel LU factorization in MUMPS ierr = MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F);CHKERRQ(ierr); ierr = MatLUFactorSymbolic(F, A, NULL, NULL, NULL);CHKERRQ(ierr); ierr = MatLUFactorNumeric(F, A, NULL);CHKERRQ(ierr); // Create spRHS Mat spRHS = NULL; // Create spRHS = spRHS^T. Two matrices that share internal matrix data structure. // Creates a new matrix object that behaves like A'. ierr = MatCreateTranspose(spRHST,&spRHS);CHKERRQ(ierr); // Get user-specified set of entries in inverse of A ierr = MatMumpsGetInverse(F,spRHS);CHKERRQ(ierr); if (displ){ ierr = PetscPrintf(PETSC_COMM_WORLD, "\n---------------\n");CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD,"First %D columns of inv(A) with sparse RHS:\n", nrhs); ierr = MatView(spRHS,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD, "---------------\n");CHKERRQ(ierr); } // Free data structures ierr = MatDestroy(&A);CHKERRQ(ierr); ierr = MatDestroy(&spRHS);CHKERRQ(ierr); ierr = MatDestroy(&spRHST);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; #endif }
