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