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