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 PetscCall(PetscInitialize(&argc,&argv,0,help)); 17 /* Create matrix and three vectors: these are all normal */ 18 PetscCall(PetscMalloc1(lda*size,&b)); 19 for (i=0; i<size; i++) { 20 for (j=0; j<size; j++) { 21 b[i+j*lda] = rand(); 22 } 23 } 24 PetscCall(MatCreate(MPI_COMM_SELF,&A)); 25 PetscCall(MatSetSizes(A,size,size,size,size)); 26 PetscCall(MatSetType(A,MATSEQDENSE)); 27 PetscCall(MatSeqDenseSetPreallocation(A,NULL)); 28 29 PetscCall(MatDenseGetArray(A,&a)); 30 for (i=0; i<size; i++) { 31 for (j=0; j<size; j++) { 32 a[i+j*size] = b[i+j*lda]; 33 } 34 } 35 PetscCall(MatDenseRestoreArray(A,&a)); 36 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 37 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 38 39 PetscCall(MatCreate(MPI_COMM_SELF,&B)); 40 PetscCall(MatSetSizes(B,size,size,size,size)); 41 PetscCall(MatSetType(B,MATSEQDENSE)); 42 PetscCall(MatSeqDenseSetPreallocation(B,b)); 43 PetscCall(MatDenseSetLDA(B,lda)); 44 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 45 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 46 47 PetscCall(PetscMalloc1(size,&x)); 48 for (i=0; i<size; i++) x[i] = 1.0; 49 PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF,1,size,x,&X)); 50 PetscCall(VecAssemblyBegin(X)); 51 PetscCall(VecAssemblyEnd(X)); 52 53 PetscCall(PetscMalloc1(size,&y)); 54 PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF,1,size,y,&Y)); 55 PetscCall(VecAssemblyBegin(Y)); 56 PetscCall(VecAssemblyEnd(Y)); 57 58 PetscCall(PetscMalloc1(size,&z)); 59 PetscCall(VecCreateSeqWithArray(MPI_COMM_SELF,1,size,z,&Z)); 60 PetscCall(VecAssemblyBegin(Z)); 61 PetscCall(VecAssemblyEnd(Z)); 62 63 /* 64 * Solve with A and B 65 */ 66 PetscCall(KSPCreate(MPI_COMM_SELF,&solver)); 67 PetscCall(KSPSetType(solver,KSPPREONLY)); 68 PetscCall(KSPGetPC(solver,&pc)); 69 PetscCall(PCSetType(pc,PCLU)); 70 PetscCall(KSPSetOperators(solver,A,A)); 71 PetscCall(KSPSolve(solver,X,Y)); 72 PetscCall(KSPSetOperators(solver,B,B)); 73 PetscCall(KSPSolve(solver,X,Z)); 74 PetscCall(VecAXPY(Z,-1.0,Y)); 75 PetscCall(VecNorm(Z,NORM_2,&nrm)); 76 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Test1; error norm=%e\n",(double)nrm)); 77 78 /* Free spaces */ 79 PetscCall(PetscFree(b)); 80 PetscCall(PetscFree(x)); 81 PetscCall(PetscFree(y)); 82 PetscCall(PetscFree(z)); 83 PetscCall(VecDestroy(&X)); 84 PetscCall(VecDestroy(&Y)); 85 PetscCall(VecDestroy(&Z)); 86 PetscCall(MatDestroy(&A)); 87 PetscCall(MatDestroy(&B)); 88 PetscCall(KSPDestroy(&solver)); 89 90 PetscCall(PetscFinalize()); 91 return 0; 92 } 93 94 /*TEST 95 96 test: 97 98 TEST*/ 99