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