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. 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);
Am Mo., 16. Dez. 2019 um 22:57 Uhr schrieb Zhang, Hong <[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 >> } >> >>
