# question

Nguyen, Hung V ERDC-ITL-MS Hung.V.Nguyen at usace.army.mil
Thu Oct 9 09:04:38 CDT 2008

```All,

I am looking for an example code that read A (in csr format) and b. Then it
builds A and b petsc format and solves Ax = b.

I found an example below, but it seems that it doesn't work.

If you have similar like an example below or let me know where is a problem,
I would appreciate very much.

Thanks,

-Hung
---
/*
Purpose: Test a sparse matrix solver.
*/
#include <stdio.h>
#include "petscksp.h"

int main(int argc,char **args)
{
/* My sample sparse matrix A */

/*
11.0   0    0  14.0  0
21.0  22.0  0  24.0  0
31.0   0   33. 34.0 35.0
0      0   43. 44.0  0
0      0    0   0   55.
*/

const int sizeMat=5;    // Matrix is 5 by 5.
int i,j;
int nonZero=12;
double val[] ={11., 14.,21., 22., 24., 31., 33., 34., 35., 43., 44.,
55.};
int col_ind[]={0, 3, 0, 1, 3, 0, 2, 3, 4, 2, 3, 4};
int row_ptr[]={0, 2, 5, 9, 11, 12};
double knB[]={2.0, 0.0, 1.0, 1.0, 2.0};
0.0363636};

// calculate row_index, vector_index and number of nonzero per row:

int nZperRow[]={3,4,2,1};
int row_ind[]={0,0, 1,1,1, 2,2,2,2, 3,3, 4};
int vec_ind[]={0,1,2,3,4};
double initX[]={9.,9.,9.,9.,9.};

/*
PetSc codes start.
*/
printf("\n*** PetSC Testing phase. ***\n");
/* Create variables of PetSc */
Vec            x,b,u;  /* approx solution, RHS, exact solution
*/ /*a
linear system, Ax = b. */
Mat            A;      /* linear system matrix */
KSP            ksp;    /* linear solver context */
PetscInt       Istart,Iend;     /* Index for local matrix of
each
processor */
PetscInt       istart,iend;     /* Index for local vector of
each
processor */
PetscViewer    viewer;
PetscMPIInt    rank;
PetscErrorCode ierr;
PetscTruth     flg;

static char help[] = "Parallel vector layout.\n\n";
/* Initialization of PetSc */
PetscInitialize(&argc,&args,(char*)0,help);
MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

/*
Create parallel matrix, specifying only its global
dimensions. :
When using MatCreate(), the matrix format can be specified at
runtime. Also, the parallel partitioning of the matrix is
determined by PETSc at runtime.

Performance tuning note:  For problems of substantial size,
preallocation of matrix memory is crucial for attaining good
performance. See the matrix chapter of the users manual for
details.
- Allocates memory for a sparse parallel matrix in AIJ format

(the default parallel PETSc format: Compressed Sparse Row).
*/
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr =
MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,sizeMat,sizeMat);CHKERRQ(ierr);
ierr = MatSetType(A, MATAIJ);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);

/*
Currently, all PETSc parallel matrix formats are partitioned
by
contiguous chunks of rows across the processors.
Determine
which
rows of the matrix are locally owned.
*/
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
printf(" Rank= %d, Istart_row= %d, Iend_row+1 = %d \n", rank, Istart,
Iend);

/*
ierr =
MatMPIAIJSetPreallocationCSR(A,row_ptr,col_ind,PETSC_NULL);CHKERRQ(ierr)
;
// Standard format, CSR
ierr =
MatSeqAIJSetPreallocation(A,0,nZperRow);CHKERRQ(ierr);
// Defining the number of nonzero for each row.
*/

ierr =
MatMPIAIJSetPreallocationCSR(A,row_ptr,col_ind,PETSC_NULL);CHKERRQ(ierr)
;
// Standard format, CSR
ierr = MatSeqAIJSetPreallocation(A,0,nZperRow);CHKERRQ(ierr); //
Defining the number of nonzero for each row.

/*
Set matrix elements in parallel.
- Each processor needs to insert only elements that it owns
locally (but any non-local elements will be sent to
the
appropriate processor during matrix assembly).
- Always specify global rows and columns of matrix entries.
*/
/* Method 1: Efficient method. */
for (i=row_ptr[Istart]; i<row_ptr[Iend]; i++)
{
//ierr =
MatSetValue(A,row_ind[i],col_ind[i],val[i],INSERT_VALUES);CHKERRQ(ierr);
ierr =
MatSetValues(A,1,&(row_ind[i]),1,&(col_ind[i]),&(val[i]),INSERT_VALUES);
CHKER
RQ(ierr);
}

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

/*
Visaualize a matrix. Set a viewer's style.
To see a dense matrix, use the following two lines:
Line1: viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD);
Line2: ierr =
PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_DENSE);
*/
ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

/*
Create parallel vectors.
*/
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,sizeMat);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
/*
PETSc parallel vectors are partitioned by
contiguous chunks of rows across the processors.
Determine
which vector are locally owned.
*/
VecGetOwnershipRange(b,&istart,&iend);
/*
Insert vector values
*/
VecSetValues(x,sizeMat,vec_ind,initX,INSERT_VALUES);
VecSetValues(b,sizeMat,vec_ind,knB,INSERT_VALUES);
/*
Assemble vector, using the 2-step process:
VecAssemblyBegin(), VecAssemblyEnd()
Computations can be done while messages are in transition
by placing code between these two statements.
*/
VecAssemblyBegin(u); VecAssemblyEnd(u);
VecAssemblyBegin(x); VecAssemblyEnd(x);
VecAssemblyBegin(b); VecAssemblyEnd(b);
/*
View the exact solution vector if desired
*/
if(rank==0) printf("Vector u: \n");
flg = 1;
if (flg) {ierr =
VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
if(rank==0) printf("Vector x: \n");
if (flg) {ierr =
VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
if(rank==0) printf("Vector b: \n");
if (flg) {ierr =
VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}

/*
Create the linear solver and set various options
*/
KSPCreate(PETSC_COMM_WORLD,&ksp);
KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
KSPSetFromOptions(ksp);

/*
Solve the linear system
*/
KSPSolve(ksp,b,x);

if(rank==0) printf("Solved Vector x: \n");
if (flg) {ierr =
VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}

/*
Free work space.  All PETSc objects should be destroyed when
they are no longer needed.
*/
ierr = KSPDestroy(ksp);CHKERRQ(ierr);
ierr = MatDestroy(A);CHKERRQ(ierr);
ierr = VecDestroy(u);CHKERRQ(ierr);  ierr =
VecDestroy(x);CHKERRQ(ierr);
ierr = VecDestroy(b);CHKERRQ(ierr);

ierr = PetscFinalize();CHKERRQ(ierr);

printf("\n");
return 0;
}

```