1 static char help[] = "Test MatMatMult(), MatTranspose(), MatTransposeMatMult() for Dense and Elemental matrices.\n\n"; 2 /* 3 Example: 4 mpiexec -n <np> ./ex104 -mat_type elemental 5 */ 6 7 #include <petscmat.h> 8 9 int main(int argc,char **argv) 10 { 11 Mat A,B,C,D; 12 PetscInt i,M=10,N=5,j,nrows,ncols,am,an,rstart,rend; 13 PetscRandom r; 14 PetscBool equal,Aiselemental; 15 PetscReal fill = 1.0; 16 IS isrows,iscols; 17 const PetscInt *rows,*cols; 18 PetscScalar *v,rval; 19 #if defined(PETSC_HAVE_ELEMENTAL) 20 PetscBool Test_MatMatMult=PETSC_TRUE; 21 #else 22 PetscBool Test_MatMatMult=PETSC_FALSE; 23 #endif 24 PetscMPIInt size; 25 26 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 27 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 28 29 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 30 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 31 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 32 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N)); 33 CHKERRQ(MatSetType(A,MATDENSE)); 34 CHKERRQ(MatSetFromOptions(A)); 35 CHKERRQ(MatSetUp(A)); 36 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&r)); 37 CHKERRQ(PetscRandomSetFromOptions(r)); 38 39 /* Set local matrix entries */ 40 CHKERRQ(MatGetOwnershipIS(A,&isrows,&iscols)); 41 CHKERRQ(ISGetLocalSize(isrows,&nrows)); 42 CHKERRQ(ISGetIndices(isrows,&rows)); 43 CHKERRQ(ISGetLocalSize(iscols,&ncols)); 44 CHKERRQ(ISGetIndices(iscols,&cols)); 45 CHKERRQ(PetscMalloc1(nrows*ncols,&v)); 46 for (i=0; i<nrows; i++) { 47 for (j=0; j<ncols; j++) { 48 CHKERRQ(PetscRandomGetValue(r,&rval)); 49 v[i*ncols+j] = rval; 50 } 51 } 52 CHKERRQ(MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES)); 53 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 54 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 55 CHKERRQ(ISRestoreIndices(isrows,&rows)); 56 CHKERRQ(ISRestoreIndices(iscols,&cols)); 57 CHKERRQ(ISDestroy(&isrows)); 58 CHKERRQ(ISDestroy(&iscols)); 59 CHKERRQ(PetscRandomDestroy(&r)); 60 61 CHKERRQ(PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&Aiselemental)); 62 63 /* Test MatCreateTranspose() and MatTranspose() */ 64 CHKERRQ(MatCreateTranspose(A,&C)); 65 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); /* B = A^T */ 66 CHKERRQ(MatMultEqual(C,B,10,&equal)); 67 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T"); 68 CHKERRQ(MatDestroy(&B)); 69 70 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&B)); 71 if (!Aiselemental) { 72 CHKERRQ(MatTranspose(B,MAT_INPLACE_MATRIX,&B)); 73 CHKERRQ(MatMultEqual(C,B,10,&equal)); 74 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"C*x != B*x"); 75 } 76 CHKERRQ(MatDestroy(&B)); 77 78 /* Test B = C*A for matrix type transpose and seqdense */ 79 if (size == 1 && !Aiselemental) { 80 CHKERRQ(MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&B)); 81 CHKERRQ(MatMatMultEqual(C,A,B,10,&equal)); 82 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"B != C*A for matrix type transpose and seqdense"); 83 CHKERRQ(MatDestroy(&B)); 84 } 85 CHKERRQ(MatDestroy(&C)); 86 87 /* Test MatMatMult() */ 88 if (Test_MatMatMult) { 89 #if !defined(PETSC_HAVE_ELEMENTAL) 90 PetscCheckFalse(size > 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"This test requires ELEMENTAL"); 91 #endif 92 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); /* B = A^T */ 93 CHKERRQ(MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C)); /* C = B*A = A^T*A */ 94 CHKERRQ(MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C)); 95 CHKERRQ(MatMatMultEqual(B,A,C,10,&equal)); 96 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x"); 97 98 /* Test MatDuplicate for matrix product */ 99 CHKERRQ(MatDuplicate(C,MAT_COPY_VALUES,&D)); 100 101 CHKERRQ(MatDestroy(&D)); 102 CHKERRQ(MatDestroy(&C)); 103 CHKERRQ(MatDestroy(&B)); 104 } 105 106 /* Test MatTransposeMatMult() */ 107 if (!Aiselemental) { 108 CHKERRQ(MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&D)); /* D = A^T*A */ 109 CHKERRQ(MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&D)); 110 CHKERRQ(MatTransposeMatMultEqual(A,A,D,10,&equal)); 111 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x"); 112 113 /* Test MatDuplicate for matrix product */ 114 CHKERRQ(MatDuplicate(D,MAT_COPY_VALUES,&C)); 115 CHKERRQ(MatDestroy(&C)); 116 CHKERRQ(MatDestroy(&D)); 117 118 /* Test D*x = A^T*C*A*x, where C is in AIJ format */ 119 CHKERRQ(MatGetLocalSize(A,&am,&an)); 120 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&C)); 121 if (size == 1) { 122 CHKERRQ(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,am,am)); 123 } else { 124 CHKERRQ(MatSetSizes(C,am,am,PETSC_DECIDE,PETSC_DECIDE)); 125 } 126 CHKERRQ(MatSetFromOptions(C)); 127 CHKERRQ(MatSetUp(C)); 128 CHKERRQ(MatGetOwnershipRange(C,&rstart,&rend)); 129 v[0] = 1.0; 130 for (i=rstart; i<rend; i++) { 131 CHKERRQ(MatSetValues(C,1,&i,1,&i,v,INSERT_VALUES)); 132 } 133 CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 134 CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 135 136 /* B = C*A, D = A^T*B */ 137 CHKERRQ(MatMatMult(C,A,MAT_INITIAL_MATRIX,1.0,&B)); 138 CHKERRQ(MatTransposeMatMult(A,B,MAT_INITIAL_MATRIX,fill,&D)); 139 CHKERRQ(MatTransposeMatMultEqual(A,B,D,10,&equal)); 140 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*B*x"); 141 142 CHKERRQ(MatDestroy(&D)); 143 CHKERRQ(MatDestroy(&C)); 144 CHKERRQ(MatDestroy(&B)); 145 } 146 147 /* Test MatMatTransposeMult() */ 148 if (!Aiselemental) { 149 PetscReal diff, scale; 150 PetscInt am, an, aM, aN; 151 152 CHKERRQ(MatGetLocalSize(A, &am, &an)); 153 CHKERRQ(MatGetSize(A, &aM, &aN)); 154 CHKERRQ(MatCreateDense(PetscObjectComm((PetscObject)A),PETSC_DECIDE, an, aM + 10, aN, NULL, &B)); 155 CHKERRQ(MatSetRandom(B, NULL)); 156 CHKERRQ(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 157 CHKERRQ(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 158 CHKERRQ(MatMatTransposeMult(A,B,MAT_INITIAL_MATRIX,fill,&D)); /* D = A*A^T */ 159 160 /* Test MatDuplicate for matrix product */ 161 CHKERRQ(MatDuplicate(D,MAT_COPY_VALUES,&C)); 162 163 CHKERRQ(MatMatTransposeMult(A,B,MAT_REUSE_MATRIX,fill,&D)); 164 CHKERRQ(MatAXPY(C, -1., D, SAME_NONZERO_PATTERN)); 165 166 CHKERRQ(MatNorm(C, NORM_FROBENIUS, &diff)); 167 CHKERRQ(MatNorm(D, NORM_FROBENIUS, &scale)); 168 PetscCheckFalse(diff > PETSC_SMALL * scale,PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatMatTransposeMult() differs between MAT_INITIAL_MATRIX and MAT_REUSE_MATRIX"); 169 CHKERRQ(MatDestroy(&C)); 170 171 CHKERRQ(MatMatTransposeMultEqual(A,B,D,10,&equal)); 172 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x"); 173 CHKERRQ(MatDestroy(&D)); 174 CHKERRQ(MatDestroy(&B)); 175 176 } 177 178 CHKERRQ(MatDestroy(&A)); 179 CHKERRQ(PetscFree(v)); 180 CHKERRQ(PetscFinalize()); 181 return 0; 182 } 183 184 /*TEST 185 186 test: 187 output_file: output/ex104.out 188 189 test: 190 suffix: 2 191 nsize: 2 192 output_file: output/ex104.out 193 194 test: 195 suffix: 3 196 nsize: 4 197 output_file: output/ex104.out 198 args: -M 23 -N 31 199 200 test: 201 suffix: 4 202 nsize: 4 203 output_file: output/ex104.out 204 args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via cyclic 205 206 test: 207 suffix: 5 208 nsize: 4 209 output_file: output/ex104.out 210 args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via allgatherv 211 212 test: 213 suffix: 6 214 args: -mat_type elemental 215 requires: elemental 216 output_file: output/ex104.out 217 218 test: 219 suffix: 7 220 nsize: 2 221 args: -mat_type elemental 222 requires: elemental 223 output_file: output/ex104.out 224 225 TEST*/ 226