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