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