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