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