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