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