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