[petsc-users] Inverse matrix empty

[Hong Zhang]

Your calling procedure is incorrect. See
~petsc/src/mat/examples/tests/ex129.c on how to use MatMatMult().
You do not need KSP object. Here are suggested calls:

  Mat F; // store matrix factors
  ierr = MatGetOrdering(Ap,MATORDERINGND,&perm,&iperm);CHKERRQ(ierr);
  ierr = MatGetFactor(Ap,MATSOLVERPETSC,MAT_FACTOR_LU,&F);CHKERRQ(ierr);
  ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr);
  info.fill = 5.0;
  ierr = MatLUFactorSymbolic(F,Ap,perm,iperm,&info);CHKERRQ(ierr);
  ierr = MatLUFactorNumeric(F,Ap,&info);CHKERRQ(ierr);

  ierr = MatMatSolve(F,Bp,Xp);CHKERRQ(ierr);

Hong

On Mon, Aug 9, 2010 at 2:14 PM, Немања Илић (Nemanja Ilic)
<nemanja.ilic.81 at gmail.com> wrote:
> Here is my code,
>
> Thank you,
>
> regards,
>
> Nemanja
>
> #include "math_petsc.h"
>
> #include <petscksp.h>
>
> #include <string.h>
>
> #include <stdio.h>
>
> #include <stdlib.h>
>
> double** getInverseMatrix(
>
> int num_rows_in,
>
> int num_columns_in,
>
> const char *file_name,
>
> const char *output_file_name)
>
> {
>
> double **X = NULL;
>
> Mat Ap; // input matrix
>
> Mat Bp, // identity matrix
>
> Xp; // result matrix
>
> PetscErrorCode ierr; // error code
>
> PetscInt i, j, // iterator variable
>
> *col, // column indices
>
> *row, // row indices
>
> ap_num_rows, // number of rows
>
> ap_num_columns; // number of columns
>
> FILE *file; // input file
>
> IS perm_rows, perm_cols;
>
> MatFactorInfo lu_info, fact_info;
>
> KSP ksp;
>
> PC pc;
>
> IS perm, iperm;
>
> MatFactorInfo info;
>
> PetscInt rank, size;
>
> ap_num_rows = num_rows_in;
>
> ap_num_columns = num_columns_in;
>
> ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
>
> ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
>
> printf("MPI_Comm_rank = %d\nMPI_Comm_size = %d\n", rank, size);
>
> // open file and read number of rows and columns
>
> ierr = PetscFOpen(PETSC_COMM_SELF, file_name, "r", &file); CHKERRQ(ierr);
>
> // alocate column and row indices
>
> col = (PetscInt*) malloc(ap_num_columns * sizeof(PetscInt));
>
> row = (PetscInt*) malloc(ap_num_rows * sizeof(PetscInt));
>
> // create matrices
>
> ierr = MatCreate(PETSC_COMM_WORLD, &Ap); CHKERRQ(ierr);
>
> ierr = MatSetSizes(Ap, PETSC_DECIDE, PETSC_DECIDE, ap_num_rows,
> ap_num_columns); CHKERRQ(ierr);
>
> ierr = MatSetFromOptions(Ap); CHKERRQ(ierr);
>
> ierr = MatCreate(PETSC_COMM_WORLD, &Bp); CHKERRQ(ierr);
>
> ierr = MatSetSizes(Bp, PETSC_DECIDE, PETSC_DECIDE, ap_num_rows,
> ap_num_columns); CHKERRQ(ierr);
>
> ierr = MatSetFromOptions(Bp); CHKERRQ(ierr);
>
> ierr = MatCreate(PETSC_COMM_WORLD, &Xp); CHKERRQ(ierr);
>
> ierr = MatSetSizes(Xp, PETSC_DECIDE, PETSC_DECIDE, ap_num_rows,
> ap_num_columns); CHKERRQ(ierr);
>
> ierr = MatSetFromOptions(Xp); CHKERRQ(ierr);
>
> /* ierr = MatCreate(PETSC_COMM_WORLD, &fact); CHKERRQ(ierr);
>
> ierr = MatSetSizes(fact, PETSC_DECIDE, PETSC_DECIDE, m, n); CHKERRQ(ierr);
>
> ierr = MatSetFromOptions(fact); CHKERRQ(ierr);
>
> */
>
> // populate row index vector
>
> for (i = 0; i < ap_num_rows; ++i)
>
> row[i] = i;
>
> // populate column index variable
>
> for (i = 0; i < ap_num_columns; ++i)
>
> col[i] = i;
>
> // now populate matrix
>
> double *row_values = (double*) malloc(ap_num_columns * sizeof(double));
>
> for (i = 0; i < ap_num_rows; ++i)
>
> {
>
> for (j = 0; j < ap_num_columns; ++j)
>
> {
>
> fscanf(file, "%lg", &row_values[j]);
>
> }
>
> printf("%f %f %f\n", row_values[j-3], row_values[j-2], row_values[j-1]);
>
> // add row to matrix
>
> ierr = MatSetValues(Ap, 1, &i, ap_num_columns, col, row_values,
> INSERT_VALUES); CHKERRQ(ierr);
>
> }
>
> ierr = MatAssemblyBegin(Ap, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> ierr = MatAssemblyEnd(Ap, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> ierr = MatView(Ap, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
>
> // factor the matrix
>
> ierr = MatGetOrdering(Ap, MATORDERING_1WD, &perm, &iperm);CHKERRQ(ierr);
>
> ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr);
>
> ierr = MatLUFactor(Ap, perm, iperm, &info);CHKERRQ(ierr);
>
> // the identity matrix
>
> double identity_values[ap_num_columns];
>
> memset(identity_values, 0, ap_num_columns * sizeof(double));
>
> for (i = 0; i < ap_num_rows; ++i)
>
> {
>
> identity_values[i] = 1;
>
> if (i != 0)
>
> identity_values[i-1] = 0;
>
> ierr = MatSetValues(Bp, 1, &i, ap_num_columns, col, identity_values,
> INSERT_VALUES); CHKERRQ(ierr);
>
> }
>
> ierr = MatAssemblyBegin(Bp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> ierr = MatAssemblyEnd(Bp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> // view identity matrix
>
> printf("This is the identity matrix:\n");
>
> ierr = MatView(Bp, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
>
> // create solver context
>
> ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); CHKERRQ(ierr);
>
> ierr = KSPSetOperators(ksp, Ap, Ap, DIFFERENT_NONZERO_PATTERN);
> CHKERRQ(ierr);
>
> // ierr = KSPSetOperators(ksp, Bp, Bp, DIFFERENT_NONZERO_PATTERN);
> CHKERRQ(ierr);
>
> // ierr = KSPSetOperators(ksp, Xp, Xp, DIFFERENT_NONZERO_PATTERN);
> CHKERRQ(ierr);
>
> ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
>
> ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
>
> ierr =
> KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
>
> ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
>
> // !!!!! get inverse matrix !!!!!
>
> ierr = MatMatSolve(Ap, Bp, Xp); CHKERRQ(ierr);
>
> // assemble matrix
>
> ierr = MatAssemblyBegin(Xp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> ierr = MatAssemblyEnd(Xp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> // print result
>
> printf("This is the result:\n");
>
> ierr = MatView(Xp, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
>
> // write result to file
>
> PetscViewer viewer;
>
> ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, output_file_name, &viewer);
> CHKERRQ(ierr);
>
> ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_COMMON);
> CHKERRQ(ierr);
>
> ierr = MatView(Xp, viewer); CHKERRQ(ierr);
>
> // clean-up
>
> ierr = MatDestroy(Ap); CHKERRQ(ierr);
>
> ierr = MatDestroy(Bp); CHKERRQ(ierr);
>
> ierr = MatDestroy(Xp); CHKERRQ(ierr);
>
> free(col);
>
> free(row);
>
> free(row_values);
>
> return X;
>
> }
>
> On Monday 09 August 2010 20:14:01 you wrote:
>
>> Send me the segment of your code for computing inv(A).
>
>> I'll check it.
>
>>
>
>> Hong
>
>>
>
>> On Mon, Aug 9, 2010 at 1:00 PM, Немања Илић (Nemanja Ilic)
>
>> <nemanja.ilic.81 at gmail.com> wrote:
>
>> > Hello,
>
>> >
>
>> > I want to get the inverse of a matrix. I use the following algorithm: I
>> > have input matrix (A), identity matrix (B) and I try to solve A * X = B,
>> > where X is the inverse matrix. I use MatMatSolve function. I included blas,
>> > fblas, and superlu libraries in my project but still I get empty inverse
>> > matrix. I factored A with LUFactor before calling MatMatSolve. What am I
>> > doing wrong?
>
>> >
>
>> > Thank you in advance,
>
>> > Regards,
>
>> > Nemanja
>
>> >
>
>>