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