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;
 15:   PC          pc;
 16:   Mat         A,B;
 17:   Vec         X,Y,Z;
 18:   MatScalar   *a;
 19:   PetscScalar *b,*x,*y,*z;
 20:   PetscReal   nrm;
 21:   PetscErrorCode ierr,size=8,lda=10, i,j;

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

 36:   MatGetArray(A,&a);
 37:   for (i=0; i<size; i++) {
 38:     for (j=0; j<size; j++) {
 39:       a[i+j*size] = b[i+j*lda];
 40:     }
 41:   }
 42:   MatRestoreArray(A,&a);
 43:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 44:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 46:   MatCreate(MPI_COMM_SELF,&B);
 47:   MatSetSizes(B,size,size,size,size);
 48:   MatSetType(B,MATSEQDENSE);
 49:   MatSeqDenseSetPreallocation(B,b);
 50:   MatSeqDenseSetLDA(B,lda);
 51:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 52:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 54:   PetscMalloc(size*sizeof(PetscScalar),&x);
 55:   for (i=0; i<size; i++) {
 56:     x[i] = 1.0;
 57:   }
 58:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,x,&X);
 59:   VecAssemblyBegin(X);
 60:   VecAssemblyEnd(X);

 62:   PetscMalloc(size*sizeof(PetscScalar),&y);
 63:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,y,&Y);
 64:   VecAssemblyBegin(Y);
 65:   VecAssemblyEnd(Y);

 67:   PetscMalloc(size*sizeof(PetscScalar),&z);
 68:   VecCreateSeqWithArray(MPI_COMM_SELF,1,size,z,&Z);
 69:   VecAssemblyBegin(Z);
 70:   VecAssemblyEnd(Z);

 72:   /*
 73:    * Solve with A and B
 74:    */
 75:   KSPCreate(MPI_COMM_SELF,&solver);
 76:   KSPSetType(solver,KSPPREONLY);
 77:   KSPGetPC(solver,&pc);
 78:   PCSetType(pc,PCLU);
 79:   KSPSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN);
 80:   KSPSolve(solver,X,Y);
 81:   KSPSetOperators(solver,B,B,DIFFERENT_NONZERO_PATTERN);
 82:   KSPSolve(solver,X,Z);
 83:   VecAXPY(Z,-1.0,Y);
 84:   VecNorm(Z,NORM_2,&nrm);
 85:   printf("Test1; error norm=%e\n",nrm);

 87:   /* Free spaces */
 88:   PetscFree(b);
 89:   PetscFree(x);
 90:   PetscFree(y);
 91:   PetscFree(z);
 92:   VecDestroy(&X);
 93:   VecDestroy(&Y);
 94:   VecDestroy(&Z);
 95:   MatDestroy(&A);
 96:   MatDestroy(&B);
 97:   KSPDestroy(&solver);

 99:   PetscFinalize();
100:   return 0;
101: }