[petsc-users] Corrrect usage of MatMumpsGetInverse()
Jan Grießer
griesser.jan at googlemail.com
Wed Dec 18 10:58:15 CST 2019
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 <hzhang at mcs.anl.gov>:
> 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
>> }
>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20191218/9903037b/attachment.html>
More information about the petsc-users
mailing list