<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; "><div><br></div> call MatView() on your Xp before doing the solve, likely that is set wrong.<div><br></div><div> Barry</div><div><br><div><div>On Aug 9, 2010, at 2:14 PM, Немања Илић (Nemanja Ilic) wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite"><div style="font-size:12pt;font-family:Sans Serif"><p>Here is my code,</p><div><br class="webkit-block-placeholder"></div><p>Thank you,</p><p>regards,</p><p>Nemanja</p><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><p>#include "math_petsc.h"</p><div><br class="webkit-block-placeholder"></div><p>#include <petscksp.h></p><div><br class="webkit-block-placeholder"></div><p>#include <string.h></p><p>#include <stdio.h></p><p>#include <stdlib.h></p><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><p> PetscInt rank, size;</p><p> ap_num_rows = num_rows_in;</p><p> ap_num_columns = num_columns_in;</p><div><br class="webkit-block-placeholder"></div><p> ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);</p><p> ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);</p><div><br class="webkit-block-placeholder"></div><p> printf("MPI_Comm_rank = %d\nMPI_Comm_size = %d\n", rank, size);</p><div><br class="webkit-block-placeholder"></div><p> // open file and read number of rows and columns</p><p> ierr = PetscFOpen(PETSC_COMM_SELF, file_name, "r", &file); CHKERRQ(ierr);</p><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><p> ierr = MatView(Ap, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);</p><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><p> // !!!!! get inverse matrix !!!!!</p><p> ierr = MatMatSolve(Ap, Bp, Xp); CHKERRQ(ierr);</p><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><p> // assemble matrix</p><p> ierr = MatAssemblyBegin(Xp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);</p><p> ierr = MatAssemblyEnd(Xp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);</p><div><br class="webkit-block-placeholder"></div><p> // print result</p><p> printf("This is the result:\n");</p><p> ierr = MatView(Xp, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);</p><div><br class="webkit-block-placeholder"></div><p> // write result to file</p><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><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><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><p> free(col);</p><p> free(row);</p><p> free(row_values);</p><div><br class="webkit-block-placeholder"></div><p> return X;</p><p>}</p><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><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>> <<a href="mailto:nemanja.ilic.81@gmail.com">nemanja.ilic.81@gmail.com</a>> 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><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div><div><br class="webkit-block-placeholder"></div>
</div></blockquote></div><br></div></body></html>