Actual source code: ex3.c

  1: /*
  2:  * Example code testing SeqDense matrices with an LDA (leading dimension
  3:  * of the user-allocated arrray) larger than M.
  4:  * This example tests the functionality of MatSolve.
  5:  */
  6: #include <stdlib.h>
 7:  #include petscmat.h
 8:  #include petscksp.h

 12: int main(int argc,char **argv)
 13: {
 14:   KSP solver; PC pc;
 15:   Mat A,B;
 16:   Vec X,Y,Z;
 17:   PetscScalar *a,*b,*x,*y,*z,one=1,mone=-1;
 18:   PetscReal nrm;
 19:   int ierr,size=8,lda=10, i,j;

 21:   PetscInitialize(&argc,&argv,0,0);

 23:   /*
 24:    * Create matrix and three vectors: these are all normal
 25:    */
 26:   PetscMalloc(size*size*sizeof(PetscScalar),&a);
 27:   PetscMalloc(lda*size*sizeof(PetscScalar),&b);
 28:   for (i=0; i<size; i++) {
 29:     for (j=0; j<size; j++) {
 30:       a[i+j*size] = rand(); b[i+j*lda] = a[i+j*size];
 31:     }
 32:   }
 33:   MatCreate(MPI_COMM_SELF,size,size,size,size,&A);
 34:   MatSetType(A,MATSEQDENSE);
 35:   MatSeqDenseSetPreallocation(A,a);
 36:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 37:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 39:   MatCreate(MPI_COMM_SELF,size,size,size,size,&B);
 40:   MatSetType(B,MATSEQDENSE);
 41:   MatSeqDenseSetPreallocation(B,b);
 42:   MatSeqDenseSetLDA(B,lda);
 43:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 44:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 46:   PetscMalloc(size*sizeof(PetscScalar),&x);
 47:   for (i=0; i<size; i++) {
 48:     x[i] = one;
 49:   }
 50:   VecCreateSeqWithArray(MPI_COMM_SELF,size,x,&X);
 51:   VecAssemblyBegin(X);
 52:   VecAssemblyEnd(X);

 54:   PetscMalloc(size*sizeof(PetscScalar),&y);
 55:   VecCreateSeqWithArray(MPI_COMM_SELF,size,y,&Y);
 56:   VecAssemblyBegin(Y);
 57:   VecAssemblyEnd(Y);

 59:   PetscMalloc(size*sizeof(PetscScalar),&z);
 60:   VecCreateSeqWithArray(MPI_COMM_SELF,size,z,&Z);
 61:   VecAssemblyBegin(Z);
 62:   VecAssemblyEnd(Z);

 64:   /*
 65:    * Solve with A and B
 66:    */
 67:   KSPCreate(MPI_COMM_SELF,&solver);
 68:   KSPSetType(solver,KSPPREONLY);
 69:   KSPGetPC(solver,&pc);
 70:   PCSetType(pc,PCLU);
 71:   KSPSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN);
 72:   KSPSetRhs(solver,X);
 73:   KSPSetSolution(solver,Y);
 74:   KSPSolve(solver);
 75:   KSPSetOperators(solver,B,B,DIFFERENT_NONZERO_PATTERN);
 76:   KSPSetSolution(solver,Z);
 77:   KSPSolve(solver);
 78:   VecAXPY(&mone,Y,Z);
 79:   VecNorm(Z,NORM_2,&nrm);
 80:   printf("Test1; error norm=%e\n",nrm);

 82:   PetscFinalize();
 83:   return 0;
 84: }