1 static char help[] = "testing SeqDense matrices with an LDA (leading dimension of the user-allocated array) larger than M.\n"; 2 /* 3 * Example code testing SeqDense matrices with an LDA (leading dimension 4 * of the user-allocated array) larger than M. 5 * This example tests the functionality of MatSetValues(), MatMult(), 6 * and MatMultTranspose(). 7 */ 8 9 #include <petscmat.h> 10 11 int main(int argc, char **argv) { 12 Mat A, A11, A12, A21, A22; 13 Vec X, X1, X2, Y, Z, Z1, Z2; 14 PetscScalar *a, *b, *x, *y, *z, v, one = 1; 15 PetscReal nrm; 16 PetscInt size = 8, size1 = 6, size2 = 2, i, j; 17 PetscRandom rnd; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &argv, 0, help)); 21 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rnd)); 22 23 /* 24 * Create matrix and three vectors: these are all normal 25 */ 26 PetscCall(PetscMalloc1(size * size, &a)); 27 PetscCall(PetscMalloc1(size * size, &b)); 28 for (i = 0; i < size; i++) { 29 for (j = 0; j < size; j++) { 30 PetscCall(PetscRandomGetValue(rnd, &a[i + j * size])); 31 b[i + j * size] = a[i + j * size]; 32 } 33 } 34 PetscCall(MatCreate(PETSC_COMM_SELF, &A)); 35 PetscCall(MatSetSizes(A, size, size, size, size)); 36 PetscCall(MatSetType(A, MATSEQDENSE)); 37 PetscCall(MatSeqDenseSetPreallocation(A, a)); 38 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 39 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 40 41 PetscCall(PetscMalloc1(size, &x)); 42 for (i = 0; i < size; i++) { x[i] = one; } 43 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, x, &X)); 44 PetscCall(VecAssemblyBegin(X)); 45 PetscCall(VecAssemblyEnd(X)); 46 47 PetscCall(PetscMalloc1(size, &y)); 48 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, y, &Y)); 49 PetscCall(VecAssemblyBegin(Y)); 50 PetscCall(VecAssemblyEnd(Y)); 51 52 PetscCall(PetscMalloc1(size, &z)); 53 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size, z, &Z)); 54 PetscCall(VecAssemblyBegin(Z)); 55 PetscCall(VecAssemblyEnd(Z)); 56 57 /* 58 * Now create submatrices and subvectors 59 */ 60 PetscCall(MatCreate(PETSC_COMM_SELF, &A11)); 61 PetscCall(MatSetSizes(A11, size1, size1, size1, size1)); 62 PetscCall(MatSetType(A11, MATSEQDENSE)); 63 PetscCall(MatSeqDenseSetPreallocation(A11, b)); 64 PetscCall(MatDenseSetLDA(A11, size)); 65 PetscCall(MatAssemblyBegin(A11, MAT_FINAL_ASSEMBLY)); 66 PetscCall(MatAssemblyEnd(A11, MAT_FINAL_ASSEMBLY)); 67 68 PetscCall(MatCreate(PETSC_COMM_SELF, &A12)); 69 PetscCall(MatSetSizes(A12, size1, size2, size1, size2)); 70 PetscCall(MatSetType(A12, MATSEQDENSE)); 71 PetscCall(MatSeqDenseSetPreallocation(A12, b + size1 * size)); 72 PetscCall(MatDenseSetLDA(A12, size)); 73 PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY)); 74 PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY)); 75 76 PetscCall(MatCreate(PETSC_COMM_SELF, &A21)); 77 PetscCall(MatSetSizes(A21, size2, size1, size2, size1)); 78 PetscCall(MatSetType(A21, MATSEQDENSE)); 79 PetscCall(MatSeqDenseSetPreallocation(A21, b + size1)); 80 PetscCall(MatDenseSetLDA(A21, size)); 81 PetscCall(MatAssemblyBegin(A21, MAT_FINAL_ASSEMBLY)); 82 PetscCall(MatAssemblyEnd(A21, MAT_FINAL_ASSEMBLY)); 83 84 PetscCall(MatCreate(PETSC_COMM_SELF, &A22)); 85 PetscCall(MatSetSizes(A22, size2, size2, size2, size2)); 86 PetscCall(MatSetType(A22, MATSEQDENSE)); 87 PetscCall(MatSeqDenseSetPreallocation(A22, b + size1 * size + size1)); 88 PetscCall(MatDenseSetLDA(A22, size)); 89 PetscCall(MatAssemblyBegin(A22, MAT_FINAL_ASSEMBLY)); 90 PetscCall(MatAssemblyEnd(A22, MAT_FINAL_ASSEMBLY)); 91 92 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size1, x, &X1)); 93 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size2, x + size1, &X2)); 94 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size1, z, &Z1)); 95 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, size2, z + size1, &Z2)); 96 97 /* 98 * Now multiple matrix times input in two ways; 99 * compare the results 100 */ 101 PetscCall(MatMult(A, X, Y)); 102 PetscCall(MatMult(A11, X1, Z1)); 103 PetscCall(MatMultAdd(A12, X2, Z1, Z1)); 104 PetscCall(MatMult(A22, X2, Z2)); 105 PetscCall(MatMultAdd(A21, X1, Z2, Z2)); 106 PetscCall(VecAXPY(Z, -1.0, Y)); 107 PetscCall(VecNorm(Z, NORM_2, &nrm)); 108 if (nrm > 100.0 * PETSC_MACHINE_EPSILON) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test1; error norm=%g\n", (double)nrm)); } 109 110 /* 111 * Next test: change both matrices 112 */ 113 PetscCall(PetscRandomGetValue(rnd, &v)); 114 i = 1; 115 j = size - 2; 116 PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES)); 117 j -= size1; 118 PetscCall(MatSetValues(A12, 1, &i, 1, &j, &v, INSERT_VALUES)); 119 PetscCall(PetscRandomGetValue(rnd, &v)); 120 i = j = size1 + 1; 121 PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES)); 122 i = j = 1; 123 PetscCall(MatSetValues(A22, 1, &i, 1, &j, &v, INSERT_VALUES)); 124 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 125 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 126 PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY)); 127 PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY)); 128 PetscCall(MatAssemblyBegin(A22, MAT_FINAL_ASSEMBLY)); 129 PetscCall(MatAssemblyEnd(A22, MAT_FINAL_ASSEMBLY)); 130 131 PetscCall(MatMult(A, X, Y)); 132 PetscCall(MatMult(A11, X1, Z1)); 133 PetscCall(MatMultAdd(A12, X2, Z1, Z1)); 134 PetscCall(MatMult(A22, X2, Z2)); 135 PetscCall(MatMultAdd(A21, X1, Z2, Z2)); 136 PetscCall(VecAXPY(Z, -1.0, Y)); 137 PetscCall(VecNorm(Z, NORM_2, &nrm)); 138 if (nrm > 100.0 * PETSC_MACHINE_EPSILON) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test2; error norm=%g\n", (double)nrm)); } 139 140 /* 141 * Transpose product 142 */ 143 PetscCall(MatMultTranspose(A, X, Y)); 144 PetscCall(MatMultTranspose(A11, X1, Z1)); 145 PetscCall(MatMultTransposeAdd(A21, X2, Z1, Z1)); 146 PetscCall(MatMultTranspose(A22, X2, Z2)); 147 PetscCall(MatMultTransposeAdd(A12, X1, Z2, Z2)); 148 PetscCall(VecAXPY(Z, -1.0, Y)); 149 PetscCall(VecNorm(Z, NORM_2, &nrm)); 150 if (nrm > 100.0 * PETSC_MACHINE_EPSILON) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test3; error norm=%g\n", (double)nrm)); } 151 PetscCall(PetscFree(a)); 152 PetscCall(PetscFree(b)); 153 PetscCall(PetscFree(x)); 154 PetscCall(PetscFree(y)); 155 PetscCall(PetscFree(z)); 156 PetscCall(MatDestroy(&A)); 157 PetscCall(MatDestroy(&A11)); 158 PetscCall(MatDestroy(&A12)); 159 PetscCall(MatDestroy(&A21)); 160 PetscCall(MatDestroy(&A22)); 161 162 PetscCall(VecDestroy(&X)); 163 PetscCall(VecDestroy(&Y)); 164 PetscCall(VecDestroy(&Z)); 165 166 PetscCall(VecDestroy(&X1)); 167 PetscCall(VecDestroy(&X2)); 168 PetscCall(VecDestroy(&Z1)); 169 PetscCall(VecDestroy(&Z2)); 170 PetscCall(PetscRandomDestroy(&rnd)); 171 PetscCall(PetscFinalize()); 172 return 0; 173 } 174 175 /*TEST 176 177 test: 178 179 TEST*/ 180