1 2 static char help[] = "Tests MatTranspose(), MatNorm(), MatAXPY() and MatAYPX().\n\n"; 3 4 #include <petscmat.h> 5 6 static PetscErrorCode TransposeAXPY(Mat C,PetscScalar alpha,Mat mat,PetscErrorCode (*f)(Mat,Mat*)) 7 { 8 Mat D,E,F,G; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 if (f == MatCreateTranspose) { 13 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: (C^T)^T = (C^T)^T + alpha * A, C=A, SAME_NONZERO_PATTERN\n");CHKERRQ(ierr); 14 } else { 15 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: (C^H)^H = (C^H)^H + alpha * A, C=A, SAME_NONZERO_PATTERN\n");CHKERRQ(ierr); 16 } 17 ierr = MatDuplicate(mat,MAT_COPY_VALUES,&C);CHKERRQ(ierr); 18 ierr = f(C,&D);CHKERRQ(ierr); 19 ierr = f(D,&E);CHKERRQ(ierr); 20 ierr = MatAXPY(E,alpha,mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 21 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 22 ierr = MatDestroy(&E);CHKERRQ(ierr); 23 ierr = MatDestroy(&D);CHKERRQ(ierr); 24 ierr = MatDestroy(&C);CHKERRQ(ierr); 25 if (f == MatCreateTranspose) { 26 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: C = C + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");CHKERRQ(ierr); 27 } else { 28 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: C = C + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");CHKERRQ(ierr); 29 } 30 ierr = MatDuplicate(mat,MAT_COPY_VALUES,&C);CHKERRQ(ierr); 31 /* MATTRANSPOSE should have a MatTranspose_Transpose or MatTranspose_HT implementation */ 32 if (f == MatCreateTranspose) { 33 ierr = MatTranspose(mat,MAT_INITIAL_MATRIX,&D);CHKERRQ(ierr); 34 } else { 35 ierr = MatHermitianTranspose(mat,MAT_INITIAL_MATRIX,&D);CHKERRQ(ierr); 36 } 37 ierr = f(D,&E);CHKERRQ(ierr); 38 ierr = MatAXPY(C,alpha,E,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 39 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 40 ierr = MatDestroy(&E);CHKERRQ(ierr); 41 ierr = MatDestroy(&D);CHKERRQ(ierr); 42 ierr = MatDestroy(&C);CHKERRQ(ierr); 43 if (f == MatCreateTranspose) { 44 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: (C^T)^T = (C^T)^T + alpha * (A^T)^T, C=A, SAME_NONZERO_PATTERN\n");CHKERRQ(ierr); 45 } else { 46 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: (C^H)^H = (C^H)^H + alpha * (A^H)^H, C=A, SAME_NONZERO_PATTERN\n");CHKERRQ(ierr); 47 } 48 ierr = MatDuplicate(mat,MAT_COPY_VALUES,&C);CHKERRQ(ierr); 49 ierr = f(C,&D);CHKERRQ(ierr); 50 ierr = f(D,&E);CHKERRQ(ierr); 51 ierr = f(mat,&F);CHKERRQ(ierr); 52 ierr = f(F,&G);CHKERRQ(ierr); 53 ierr = MatAXPY(E,alpha,G,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 54 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 55 ierr = MatDestroy(&G);CHKERRQ(ierr); 56 ierr = MatDestroy(&F);CHKERRQ(ierr); 57 ierr = MatDestroy(&E);CHKERRQ(ierr); 58 ierr = MatDestroy(&D);CHKERRQ(ierr); 59 ierr = MatDestroy(&C);CHKERRQ(ierr); 60 PetscFunctionReturn(0); 61 } 62 63 int main(int argc,char **argv) 64 { 65 Mat mat,tmat = 0; 66 PetscInt m = 7,n,i,j,rstart,rend,rect = 0; 67 PetscErrorCode ierr; 68 PetscMPIInt size,rank; 69 PetscBool flg; 70 PetscScalar v, alpha; 71 PetscReal normf,normi,norm1; 72 73 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 74 ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);CHKERRQ(ierr); 75 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 76 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 77 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 78 n = m; 79 ierr = PetscOptionsHasName(NULL,NULL,"-rectA",&flg);CHKERRQ(ierr); 80 if (flg) {n += 2; rect = 1;} 81 ierr = PetscOptionsHasName(NULL,NULL,"-rectB",&flg);CHKERRQ(ierr); 82 if (flg) {n -= 2; rect = 1;} 83 84 /* ------- Assemble matrix, test MatValid() --------- */ 85 ierr = MatCreate(PETSC_COMM_WORLD,&mat);CHKERRQ(ierr); 86 ierr = MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 87 ierr = MatSetFromOptions(mat);CHKERRQ(ierr); 88 ierr = MatSetUp(mat);CHKERRQ(ierr); 89 ierr = MatGetOwnershipRange(mat,&rstart,&rend);CHKERRQ(ierr); 90 for (i=rstart; i<rend; i++) { 91 for (j=0; j<n; j++) { 92 v = 10.0*i+j; 93 ierr = MatSetValues(mat,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 94 } 95 } 96 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 97 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 98 99 /* ----------------- Test MatNorm() ----------------- */ 100 ierr = MatNorm(mat,NORM_FROBENIUS,&normf);CHKERRQ(ierr); 101 ierr = MatNorm(mat,NORM_1,&norm1);CHKERRQ(ierr); 102 ierr = MatNorm(mat,NORM_INFINITY,&normi);CHKERRQ(ierr); 103 ierr = PetscPrintf(PETSC_COMM_WORLD,"original A: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi);CHKERRQ(ierr); 104 ierr = MatView(mat,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 105 106 /* --------------- Test MatTranspose() -------------- */ 107 ierr = PetscOptionsHasName(NULL,NULL,"-in_place",&flg);CHKERRQ(ierr); 108 if (!rect && flg) { 109 ierr = MatTranspose(mat,MAT_REUSE_MATRIX,&mat);CHKERRQ(ierr); /* in-place transpose */ 110 tmat = mat; mat = 0; 111 } else { /* out-of-place transpose */ 112 ierr = MatTranspose(mat,MAT_INITIAL_MATRIX,&tmat);CHKERRQ(ierr); 113 } 114 115 /* ----------------- Test MatNorm() ----------------- */ 116 /* Print info about transpose matrix */ 117 ierr = MatNorm(tmat,NORM_FROBENIUS,&normf);CHKERRQ(ierr); 118 ierr = MatNorm(tmat,NORM_1,&norm1);CHKERRQ(ierr); 119 ierr = MatNorm(tmat,NORM_INFINITY,&normi);CHKERRQ(ierr); 120 ierr = PetscPrintf(PETSC_COMM_WORLD,"B = A^T: Frobenious norm = %g, one norm = %g, infinity norm = %g\n",(double)normf,(double)norm1,(double)normi);CHKERRQ(ierr); 121 ierr = MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 122 123 /* ----------------- Test MatAXPY(), MatAYPX() ----------------- */ 124 if (mat && !rect) { 125 alpha = 1.0; 126 ierr = PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL);CHKERRQ(ierr); 127 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: B = B + alpha * A\n");CHKERRQ(ierr); 128 ierr = MatAXPY(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 129 ierr = MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 130 131 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAYPX: B = alpha*B + A\n");CHKERRQ(ierr); 132 ierr = MatAYPX(tmat,alpha,mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 133 ierr = MatView(tmat,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 134 } 135 136 { 137 Mat C; 138 alpha = 1.0; 139 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: C = C + alpha * A, C=A, SAME_NONZERO_PATTERN\n");CHKERRQ(ierr); 140 ierr = MatDuplicate(mat,MAT_COPY_VALUES,&C);CHKERRQ(ierr); 141 ierr = MatAXPY(C,alpha,mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 142 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 143 ierr = MatDestroy(&C);CHKERRQ(ierr); 144 ierr = TransposeAXPY(C,alpha,mat,MatCreateTranspose);CHKERRQ(ierr); 145 ierr = TransposeAXPY(C,alpha,mat,MatCreateHermitianTranspose);CHKERRQ(ierr); 146 } 147 148 { 149 Mat matB; 150 /* get matB that has nonzeros of mat in all even numbers of row and col */ 151 ierr = MatCreate(PETSC_COMM_WORLD,&matB);CHKERRQ(ierr); 152 ierr = MatSetSizes(matB,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 153 ierr = MatSetFromOptions(matB);CHKERRQ(ierr); 154 ierr = MatSetUp(matB);CHKERRQ(ierr); 155 ierr = MatGetOwnershipRange(matB,&rstart,&rend);CHKERRQ(ierr); 156 if (rstart % 2 != 0) rstart++; 157 for (i=rstart; i<rend; i += 2) { 158 for (j=0; j<n; j += 2) { 159 v = 10.0*i+j; 160 ierr = MatSetValues(matB,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 161 } 162 } 163 ierr = MatAssemblyBegin(matB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 164 ierr = MatAssemblyEnd(matB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 165 PetscPrintf(PETSC_COMM_WORLD," A: original matrix:\n"); 166 ierr = MatView(mat,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 167 PetscPrintf(PETSC_COMM_WORLD," B(a subset of A):\n"); 168 ierr = MatView(matB,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 169 ierr = PetscPrintf(PETSC_COMM_WORLD,"MatAXPY: B = B + alpha * A, SUBSET_NONZERO_PATTERN\n");CHKERRQ(ierr); 170 ierr = MatAXPY(mat,alpha,matB,SUBSET_NONZERO_PATTERN);CHKERRQ(ierr); 171 ierr = MatView(mat,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 172 ierr = MatDestroy(&matB);CHKERRQ(ierr); 173 } 174 175 ierr = PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 176 /* Free data structures */ 177 ierr = MatDestroy(&mat);CHKERRQ(ierr); 178 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 179 ierr = PetscFinalize(); 180 return ierr; 181 } 182 183 184 185 186 187 /*TEST 188 189 test: 190 suffix: 11_A 191 args: -mat_type seqaij -rectA 192 filter: grep -v "Mat Object" 193 194 test: 195 suffix: 12_A 196 args: -mat_type seqdense -rectA 197 filter: grep -v "Mat Object" 198 199 test: 200 suffix: 11_B 201 args: -mat_type seqaij -rectB 202 filter: grep -v "Mat Object" 203 204 test: 205 suffix: 12_B 206 args: -mat_type seqdense -rectB 207 filter: grep -v "Mat Object" 208 209 test: 210 suffix: 21 211 args: -mat_type mpiaij 212 213 test: 214 suffix: 22 215 args: -mat_type mpidense 216 217 test: 218 suffix: 23 219 nsize: 3 220 args: -mat_type mpiaij 221 filter: grep -v type | grep -v "MPI processes" 222 223 test: 224 suffix: 24 225 nsize: 3 226 args: -mat_type mpidense 227 228 test: 229 suffix: 2_aijcusparse_1 230 args: -mat_type mpiaijcusparse 231 output_file: output/ex2_23.out 232 requires: cuda 233 filter: grep -v type | grep -v "MPI processes" 234 235 test: 236 suffix: 2_aijcusparse_2 237 nsize: 3 238 args: -mat_type mpiaijcusparse 239 output_file: output/ex2_23.out 240 requires: cuda 241 filter: grep -v type | grep -v "MPI processes" 242 243 test: 244 suffix: 3 245 nsize: 2 246 args: -mat_type mpiaij -rectA 247 248 test: 249 suffix: 3_aijcusparse 250 nsize: 2 251 args: -mat_type mpiaijcusparse -rectA 252 requires: cuda 253 254 test: 255 suffix: 4 256 nsize: 2 257 args: -mat_type mpidense -rectA 258 259 test: 260 suffix: aijcusparse_1 261 args: -mat_type seqaijcusparse -rectA 262 filter: grep -v "Mat Object" 263 output_file: output/ex2_11_A_aijcusparse.out 264 requires: cuda 265 266 test: 267 suffix: aijcusparse_2 268 args: -mat_type seqaijcusparse -rectB 269 filter: grep -v "Mat Object" 270 output_file: output/ex2_11_B_aijcusparse.out 271 requires: cuda 272 273 TEST*/ 274