[petsc-users] Inverse matrix empty

[Barry Smith]

  call MatView() on your Xp before doing the solve, likely that is set wrong.

  Barry

On Aug 9, 2010, at 2:14 PM, Немања Илић (Nemanja Ilic) 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
> 
> > >
> 
> >
> 
> 
> 
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20100809/7ad4f064/attachment-0001.htm>