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