1 static char help[] = " * Example code testing SeqDense matrices with an LDA (leading dimension of the user-allocated array) larger than M.\n"; 2 3 #include <petscksp.h> 4 5 int main(int argc,char **argv) 6 { 7 KSP solver; 8 PC pc; 9 Mat A,B; 10 Vec X,Y,Z; 11 MatScalar *a; 12 PetscScalar *b,*x,*y,*z; 13 PetscReal nrm; 14 PetscInt size = 8,lda=10, i,j; 15 16 PetscFunctionBeginUser; 17 PetscCall(PetscInitialize(&argc,&argv,0,help)); 18 /* Create matrix and three vectors: these are all normal */ 19 PetscCall(PetscMalloc1(lda*size,&b)); 20 for (i=0; i<size; i++) { 21 for (j=0; j<size; j++) { 22 b[i+j*lda] = rand(); 23 } 24 } 25 PetscCall(MatCreate(MPI_COMM_SELF,&A)); 26 PetscCall(MatSetSizes(A,size,size,size,size)); 27 PetscCall(MatSetType(A,MATSEQDENSE)); 28 PetscCall(MatSeqDenseSetPreallocation(A,NULL)); 29 30 PetscCall(MatDenseGetArray(A,&a)); 31 for (i=0; i<size; i++) { 32 for (j=0; j<size; j++) { 33 a[i+j*size] = b[i+j*lda]; 34 } 35 } 36 PetscCall(MatDenseRestoreArray(A,&a)); 37 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 38 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 39 40 PetscCall(MatCreate(MPI_COMM_SELF,&B)); 41 PetscCall(MatSetSizes(B,size,size,size,size)); 42 PetscCall(MatSetType(B,MATSEQDENSE)); 43 PetscCall(MatSeqDenseSetPreallocation(B,b)); 44 PetscCall(MatDenseSetLDA(B,lda)); 45 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 46 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 47 48 PetscCall(PetscMalloc1(size,&x)); 49 for (i=0; i<size; i++) x[i] = 1.0; 50 PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF,1,size,x,&X)); 51 PetscCall(VecAssemblyBegin(X)); 52 PetscCall(VecAssemblyEnd(X)); 53 54 PetscCall(PetscMalloc1(size,&y)); 55 PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF,1,size,y,&Y)); 56 PetscCall(VecAssemblyBegin(Y)); 57 PetscCall(VecAssemblyEnd(Y)); 58 59 PetscCall(PetscMalloc1(size,&z)); 60 PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF,1,size,z,&Z)); 61 PetscCall(VecAssemblyBegin(Z)); 62 PetscCall(VecAssemblyEnd(Z)); 63 64 /* 65 * Solve with A and B 66 */ 67 PetscCall(KSPCreate(MPI_COMM_SELF,&solver)); 68 PetscCall(KSPSetType(solver,KSPPREONLY)); 69 PetscCall(KSPGetPC(solver,&pc)); 70 PetscCall(PCSetType(pc,PCLU)); 71 PetscCall(KSPSetOperators(solver,A,A)); 72 PetscCall(KSPSolve(solver,X,Y)); 73 PetscCall(KSPSetOperators(solver,B,B)); 74 PetscCall(KSPSolve(solver,X,Z)); 75 PetscCall(VecAXPY(Z,-1.0,Y)); 76 PetscCall(VecNorm(Z,NORM_2,&nrm)); 77 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Test1; error norm=%e\n",(double)nrm)); 78 79 /* Free spaces */ 80 PetscCall(PetscFree(b)); 81 PetscCall(PetscFree(x)); 82 PetscCall(PetscFree(y)); 83 PetscCall(PetscFree(z)); 84 PetscCall(VecDestroy(&X)); 85 PetscCall(VecDestroy(&Y)); 86 PetscCall(VecDestroy(&Z)); 87 PetscCall(MatDestroy(&A)); 88 PetscCall(MatDestroy(&B)); 89 PetscCall(KSPDestroy(&solver)); 90 91 PetscCall(PetscFinalize()); 92 return 0; 93 } 94 95 /*TEST 96 97 test: 98 99 TEST*/ 100