1 2 static char help[] = "Tests MatGetRowMax(), MatGetRowMin(), MatGetRowMaxAbs()\n"; 3 4 #include <petscmat.h> 5 6 #define M 5 7 #define N 6 8 9 int main(int argc,char **args) 10 { 11 Mat A; 12 Vec min,max,maxabs,e; 13 PetscInt m,n,j,imin[M],imax[M],imaxabs[M],indices[N],row,testcase=0; 14 PetscScalar values[N]; 15 PetscMPIInt size,rank; 16 PetscReal enorm; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 20 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 21 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 22 PetscCall(PetscOptionsGetInt(NULL,NULL,"-testcase",&testcase,NULL)); 23 24 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 25 if (testcase == 1) { /* proc[0] holds entire A and other processes have no entry */ 26 if (rank == 0) { 27 PetscCall(MatSetSizes(A,M,N,PETSC_DECIDE,PETSC_DECIDE)); 28 } else { 29 PetscCall(MatSetSizes(A,0,0,PETSC_DECIDE,PETSC_DECIDE)); 30 } 31 testcase = 0; 32 } else { 33 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N)); 34 } 35 PetscCall(MatSetFromOptions(A)); 36 PetscCall(MatSetUp(A)); 37 38 if (rank == 0) { /* proc[0] sets matrix A */ 39 for (j=0; j<N; j++) indices[j] = j; 40 switch (testcase) { 41 case 1: /* see testcast 0 */ 42 break; 43 case 2: 44 row = 0; 45 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -4.0; values[4] = 1.0; values[5] = 1.0; 46 PetscCall(MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES)); 47 row = 2; 48 indices[0] = 0; indices[1] = 3; indices[2] = 5; 49 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 50 PetscCall(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 51 row = 3; 52 indices[0] = 0; indices[1] = 1; indices[2] = 4; 53 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 54 PetscCall(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 55 row = 4; 56 indices[0] = 0; indices[1] = 1; indices[2] = 2; 57 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 58 PetscCall(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 59 break; 60 case 3: 61 row = 0; 62 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 63 PetscCall(MatSetValues(A,1,&row,3,indices+1,values,INSERT_VALUES)); 64 row = 1; 65 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 66 PetscCall(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 67 row = 2; 68 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; 69 PetscCall(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 70 row = 3; 71 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -1.0; 72 PetscCall(MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES)); 73 row = 4; 74 values[0] = -2.0; values[1] = -2.0; values[2] = -2.0; values[3] = -1.0; 75 PetscCall(MatSetValues(A,1,&row,4,indices,values,INSERT_VALUES)); 76 break; 77 78 default: 79 row = 0; 80 values[0] = -1.0; values[1] = 0.0; values[2] = 1.0; values[3] = 3.0; values[4] = 4.0; values[5] = -5.0; 81 PetscCall(MatSetValues(A,1,&row,N,indices,values,INSERT_VALUES)); 82 row = 1; 83 PetscCall(MatSetValues(A,1,&row,3,indices,values,INSERT_VALUES)); 84 row = 3; 85 PetscCall(MatSetValues(A,1,&row,1,indices+4,values+4,INSERT_VALUES)); 86 row = 4; 87 PetscCall(MatSetValues(A,1,&row,2,indices+4,values+4,INSERT_VALUES)); 88 } 89 } 90 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 91 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 92 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 93 94 PetscCall(MatGetLocalSize(A, &m,&n)); 95 PetscCall(VecCreate(PETSC_COMM_WORLD,&min)); 96 PetscCall(VecSetSizes(min,m,PETSC_DECIDE)); 97 PetscCall(VecSetFromOptions(min)); 98 PetscCall(VecDuplicate(min,&max)); 99 PetscCall(VecDuplicate(min,&maxabs)); 100 PetscCall(VecDuplicate(min,&e)); 101 102 /* MatGetRowMax() */ 103 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMax\n")); 104 PetscCall(MatGetRowMax(A,max,NULL)); 105 PetscCall(MatGetRowMax(A,max,imax)); 106 PetscCall(VecView(max,PETSC_VIEWER_STDOUT_WORLD)); 107 PetscCall(VecGetLocalSize(max,&n)); 108 PetscCall(PetscIntView(n,imax,PETSC_VIEWER_STDOUT_WORLD)); 109 110 /* MatGetRowMin() */ 111 PetscCall(MatScale(A,-1.0)); 112 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 113 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMin\n")); 114 PetscCall(MatGetRowMin(A,min,NULL)); 115 PetscCall(MatGetRowMin(A,min,imin)); 116 117 PetscCall(VecWAXPY(e,1.0,max,min)); /* e = max + min */ 118 PetscCall(VecNorm(e,NORM_INFINITY,&enorm)); 119 PetscCheck(enorm <= PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"max+min > PETSC_MACHINE_EPSILON "); 120 for (j = 0; j < n; j++) { 121 PetscCheck(imin[j] == imax[j],PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT,j,imin[j],imax[j]); 122 } 123 124 /* MatGetRowMaxAbs() */ 125 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs\n")); 126 PetscCall(MatGetRowMaxAbs(A,maxabs,NULL)); 127 PetscCall(MatGetRowMaxAbs(A,maxabs,imaxabs)); 128 PetscCall(VecView(maxabs,PETSC_VIEWER_STDOUT_WORLD)); 129 PetscCall(PetscIntView(n,imaxabs,PETSC_VIEWER_STDOUT_WORLD)); 130 131 /* MatGetRowMinAbs() */ 132 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMinAbs\n")); 133 PetscCall(MatGetRowMinAbs(A,min,NULL)); 134 PetscCall(MatGetRowMinAbs(A,min,imin)); 135 PetscCall(VecView(min,PETSC_VIEWER_STDOUT_WORLD)); 136 PetscCall(PetscIntView(n,imin,PETSC_VIEWER_STDOUT_WORLD)); 137 138 if (size == 1) { 139 /* Test MatGetRowMax, MatGetRowMin and MatGetRowMaxAbs for SeqDense and MPIBAIJ matrix */ 140 Mat Adense; 141 Vec max_d,maxabs_d; 142 PetscCall(VecDuplicate(min,&max_d)); 143 PetscCall(VecDuplicate(min,&maxabs_d)); 144 145 PetscCall(MatScale(A,-1.0)); 146 PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense)); 147 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMax for seqdense matrix\n")); 148 PetscCall(MatGetRowMax(Adense,max_d,imax)); 149 150 PetscCall(VecWAXPY(e,-1.0,max,max_d)); /* e = -max + max_d */ 151 PetscCall(VecNorm(e,NORM_INFINITY,&enorm)); 152 PetscCheck(enorm <= PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-max + max_d) %g > PETSC_MACHINE_EPSILON",(double)enorm); 153 154 PetscCall(MatScale(Adense,-1.0)); 155 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMin for seqdense matrix\n")); 156 PetscCall(MatGetRowMin(Adense,min,imin)); 157 158 PetscCall(VecWAXPY(e,1.0,max,min)); /* e = max + min */ 159 PetscCall(VecNorm(e,NORM_INFINITY,&enorm)); 160 PetscCheck(enorm <= PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"max+min > PETSC_MACHINE_EPSILON "); 161 for (j = 0; j < n; j++) { 162 PetscCheck(imin[j] == imax[j],PETSC_COMM_SELF,PETSC_ERR_PLIB,"imin[%" PetscInt_FMT "] %" PetscInt_FMT " != imax %" PetscInt_FMT " for seqdense matrix",j,imin[j],imax[j]); 163 } 164 165 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"MatGetRowMaxAbs for seqdense matrix\n")); 166 PetscCall(MatGetRowMaxAbs(Adense,maxabs_d,imaxabs)); 167 PetscCall(VecWAXPY(e,-1.0,maxabs,maxabs_d)); /* e = -maxabs + maxabs_d */ 168 PetscCall(VecNorm(e,NORM_INFINITY,&enorm)); 169 PetscCheck(enorm <= PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",(double)enorm); 170 171 PetscCall(MatDestroy(&Adense)); 172 PetscCall(VecDestroy(&max_d)); 173 PetscCall(VecDestroy(&maxabs_d)); 174 } 175 176 { /* BAIJ matrix */ 177 Mat B; 178 Vec maxabsB,maxabsB2; 179 PetscInt bs=2,*imaxabsB,*imaxabsB2,rstart,rend,cstart,cend,ncols,col,Brows[2],Bcols[2]; 180 const PetscInt *cols; 181 const PetscScalar *vals,*vals2; 182 PetscScalar Bvals[4]; 183 184 PetscCall(PetscMalloc2(M,&imaxabsB,bs*M,&imaxabsB2)); 185 186 /* bs = 1 */ 187 PetscCall(MatConvert(A,MATMPIBAIJ,MAT_INITIAL_MATRIX,&B)); 188 PetscCall(VecDuplicate(min,&maxabsB)); 189 PetscCall(MatGetRowMaxAbs(B,maxabsB,NULL)); 190 PetscCall(MatGetRowMaxAbs(B,maxabsB,imaxabsB)); 191 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n MatGetRowMaxAbs for BAIJ matrix\n")); 192 PetscCall(VecWAXPY(e,-1.0,maxabs,maxabsB)); /* e = -maxabs + maxabsB */ 193 PetscCall(VecNorm(e,NORM_INFINITY,&enorm)); 194 PetscCheck(enorm <= PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"norm(-maxabs + maxabs_d) %g > PETSC_MACHINE_EPSILON",(double)enorm); 195 196 for (j = 0; j < n; j++) { 197 PetscCheck(imaxabs[j] == imaxabsB[j],PETSC_COMM_SELF,PETSC_ERR_PLIB,"imaxabs[%" PetscInt_FMT "] %" PetscInt_FMT " != imaxabsB %" PetscInt_FMT,j,imin[j],imax[j]); 198 } 199 PetscCall(MatDestroy(&B)); 200 201 /* Test bs = 2: Create B with bs*bs block structure of A */ 202 PetscCall(VecCreate(PETSC_COMM_WORLD,&maxabsB2)); 203 PetscCall(VecSetSizes(maxabsB2,bs*m,PETSC_DECIDE)); 204 PetscCall(VecSetFromOptions(maxabsB2)); 205 206 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 207 PetscCall(MatGetOwnershipRangeColumn(A,&cstart,&cend)); 208 PetscCall(MatCreate(PETSC_COMM_WORLD,&B)); 209 PetscCall(MatSetSizes(B,bs*(rend-rstart),bs*(cend-cstart),PETSC_DECIDE,PETSC_DECIDE)); 210 PetscCall(MatSetFromOptions(B)); 211 PetscCall(MatSetUp(B)); 212 213 for (row=rstart; row<rend; row++) { 214 PetscCall(MatGetRow(A,row,&ncols,&cols,&vals)); 215 for (col=0; col<ncols; col++) { 216 for (j=0; j<bs; j++) { 217 Brows[j] = bs*row + j; 218 Bcols[j] = bs*cols[col]+j; 219 } 220 for (j=0; j<bs*bs; j++) Bvals[j] = vals[col]; 221 PetscCall(MatSetValues(B,bs,Brows,bs,Bcols,Bvals,INSERT_VALUES)); 222 } 223 PetscCall(MatRestoreRow(A,row,&ncols,&cols,&vals)); 224 } 225 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 226 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 227 228 PetscCall(MatGetRowMaxAbs(B,maxabsB2,imaxabsB2)); 229 230 /* Check maxabsB2 and imaxabsB2 */ 231 PetscCall(VecGetArrayRead(maxabsB,&vals)); 232 PetscCall(VecGetArrayRead(maxabsB2,&vals2)); 233 for (row=0; row<m; row++) { 234 PetscCheck(PetscAbsScalar(vals[row] - vals2[bs*row]) <= PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_PLIB,"row %" PetscInt_FMT " maxabsB != maxabsB2",row); 235 } 236 PetscCall(VecRestoreArrayRead(maxabsB,&vals)); 237 PetscCall(VecRestoreArrayRead(maxabsB2,&vals2)); 238 239 for (col=0; col<n; col++) { 240 PetscCheck(imaxabsB[col] == imaxabsB2[bs*col]/bs,PETSC_COMM_SELF,PETSC_ERR_PLIB,"col %" PetscInt_FMT " imaxabsB != imaxabsB2",col); 241 } 242 PetscCall(VecDestroy(&maxabsB)); 243 PetscCall(MatDestroy(&B)); 244 PetscCall(VecDestroy(&maxabsB2)); 245 PetscCall(PetscFree2(imaxabsB,imaxabsB2)); 246 } 247 248 PetscCall(VecDestroy(&min)); 249 PetscCall(VecDestroy(&max)); 250 PetscCall(VecDestroy(&maxabs)); 251 PetscCall(VecDestroy(&e)); 252 PetscCall(MatDestroy(&A)); 253 PetscCall(PetscFinalize()); 254 return 0; 255 } 256 257 /*TEST 258 259 test: 260 output_file: output/ex114.out 261 262 test: 263 suffix: 2 264 args: -testcase 1 265 output_file: output/ex114.out 266 267 test: 268 suffix: 3 269 args: -testcase 2 270 output_file: output/ex114_3.out 271 272 test: 273 suffix: 4 274 args: -testcase 3 275 output_file: output/ex114_4.out 276 277 test: 278 suffix: 5 279 nsize: 3 280 args: -testcase 0 281 output_file: output/ex114_5.out 282 283 test: 284 suffix: 6 285 nsize: 3 286 args: -testcase 1 287 output_file: output/ex114_6.out 288 289 test: 290 suffix: 7 291 nsize: 3 292 args: -testcase 2 293 output_file: output/ex114_7.out 294 295 test: 296 suffix: 8 297 nsize: 3 298 args: -testcase 3 299 output_file: output/ex114_8.out 300 301 TEST*/ 302