1 2 static char help[] = "Tests various routines in MatMPISBAIJ format.\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Vec x,y,u,s1,s2; 9 Mat A,sA,sB; 10 PetscRandom rctx; 11 PetscReal r1,r2,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON; 12 PetscScalar one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1; 13 PetscInt n,col[3],n1,block,row,i,j,i2,j2,Ii,J,rstart,rend,bs=1,mbs=16,d_nz=3,o_nz=3,prob=2; 14 PetscErrorCode ierr; 15 PetscMPIInt size,rank; 16 PetscBool flg; 17 MatType type; 18 19 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 20 ierr = PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);CHKERRQ(ierr); 21 ierr = PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);CHKERRQ(ierr); 22 23 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 24 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 25 26 n = mbs*bs; 27 28 /* Assemble MPISBAIJ matrix sA */ 29 ierr = MatCreate(PETSC_COMM_WORLD,&sA);CHKERRQ(ierr); 30 ierr = MatSetSizes(sA,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr); 31 ierr = MatSetType(sA,MATSBAIJ);CHKERRQ(ierr); 32 ierr = MatSetFromOptions(sA);CHKERRQ(ierr); 33 ierr = MatGetType(sA,&type);CHKERRQ(ierr); 34 ierr = MatMPISBAIJSetPreallocation(sA,bs,d_nz,NULL,o_nz,NULL);CHKERRQ(ierr); 35 ierr = MatSeqSBAIJSetPreallocation(sA,bs,d_nz,NULL);CHKERRQ(ierr); 36 ierr = MatSetOption(sA,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 37 38 if (bs == 1) { 39 if (prob == 1) { /* tridiagonal matrix */ 40 value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; 41 for (i=1; i<n-1; i++) { 42 col[0] = i-1; col[1] = i; col[2] = i+1; 43 ierr = MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 44 } 45 i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1; 46 value[0]= 0.1; value[1]=-1; value[2]=2; 47 ierr = MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 48 49 i = 0; col[0] = 0; col[1] = 1; col[2]=n-1; 50 value[0] = 2.0; value[1] = -1.0; value[2]=0.1; 51 ierr = MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 52 } else if (prob ==2) { /* matrix for the five point stencil */ 53 n1 = (int) PetscSqrtReal((PetscReal)n); 54 if (n1*n1 != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1"); 55 56 for (i=0; i<n1; i++) { 57 for (j=0; j<n1; j++) { 58 Ii = j + n1*i; 59 if (i>0) {J = Ii - n1; ierr = MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);} 60 if (i<n1-1) {J = Ii + n1; ierr = MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);} 61 if (j>0) {J = Ii - 1; ierr = MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);} 62 if (j<n1-1) {J = Ii + 1; ierr = MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);} 63 ierr = MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);CHKERRQ(ierr); 64 } 65 } 66 } 67 /* end of if (bs == 1) */ 68 } else { /* bs > 1 */ 69 for (block=0; block<n/bs; block++) { 70 /* diagonal blocks */ 71 value[0] = -1.0; value[1] = 4.0; value[2] = -1.0; 72 for (i=1+block*bs; i<bs-1+block*bs; i++) { 73 col[0] = i-1; col[1] = i; col[2] = i+1; 74 ierr = MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 75 } 76 i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs; 77 value[0]=-1.0; value[1]=4.0; 78 ierr = MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 79 80 i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs; 81 value[0]=4.0; value[1] = -1.0; 82 ierr = MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 83 } 84 /* off-diagonal blocks */ 85 value[0]=-1.0; 86 for (i=0; i<(n/bs-1)*bs; i++) { 87 col[0]=i+bs; 88 ierr = MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 89 col[0]=i; row=i+bs; 90 ierr = MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 91 } 92 } 93 ierr = MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 94 ierr = MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 95 96 /* Test MatView() */ 97 ierr = MatCreateBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,NULL,o_nz,NULL,&A);CHKERRQ(ierr); 98 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 99 100 if (bs == 1) { 101 if (prob == 1) { /* tridiagonal matrix */ 102 value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; 103 for (i=1; i<n-1; i++) { 104 col[0] = i-1; col[1] = i; col[2] = i+1; 105 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 106 } 107 i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1; 108 value[0]= 0.1; value[1]=-1; value[2]=2; 109 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 110 111 i = 0; col[0] = 0; col[1] = 1; col[2]=n-1; 112 value[0] = 2.0; value[1] = -1.0; value[2]=0.1; 113 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 114 } else if (prob ==2) { /* matrix for the five point stencil */ 115 n1 = (int) PetscSqrtReal((PetscReal)n); 116 for (i=0; i<n1; i++) { 117 for (j=0; j<n1; j++) { 118 Ii = j + n1*i; 119 if (i>0) {J = Ii - n1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);} 120 if (i<n1-1) {J = Ii + n1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);} 121 if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);} 122 if (j<n1-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);} 123 ierr = MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);CHKERRQ(ierr); 124 } 125 } 126 } 127 /* end of if (bs == 1) */ 128 } else { /* bs > 1 */ 129 for (block=0; block<n/bs; block++) { 130 /* diagonal blocks */ 131 value[0] = -1.0; value[1] = 4.0; value[2] = -1.0; 132 for (i=1+block*bs; i<bs-1+block*bs; i++) { 133 col[0] = i-1; col[1] = i; col[2] = i+1; 134 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 135 } 136 i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs; 137 value[0]=-1.0; value[1]=4.0; 138 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 139 140 i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs; 141 value[0]=4.0; value[1] = -1.0; 142 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 143 } 144 /* off-diagonal blocks */ 145 value[0]=-1.0; 146 for (i=0; i<(n/bs-1)*bs; i++) { 147 col[0]=i+bs; 148 ierr = MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 149 col[0]=i; row=i+bs; 150 ierr = MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 151 } 152 } 153 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 154 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 155 156 /* Test MatGetSize(), MatGetLocalSize() */ 157 ierr = MatGetSize(sA, &i,&j); CHKERRQ(ierr); 158 ierr = MatGetSize(A, &i2,&j2);CHKERRQ(ierr); 159 i -= i2; j -= j2; 160 if (i || j) { 161 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);CHKERRQ(ierr); 162 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 163 } 164 165 ierr = MatGetLocalSize(sA, &i,&j);CHKERRQ(ierr); 166 ierr = MatGetLocalSize(A, &i2,&j2);CHKERRQ(ierr); 167 i2 -= i; j2 -= j; 168 if (i2 || j2) { 169 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);CHKERRQ(ierr); 170 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 171 } 172 173 /* vectors */ 174 /*--------------------*/ 175 /* i is obtained from MatGetLocalSize() */ 176 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 177 ierr = VecSetSizes(x,i,PETSC_DECIDE);CHKERRQ(ierr); 178 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 179 ierr = VecDuplicate(x,&y);CHKERRQ(ierr); 180 ierr = VecDuplicate(x,&u);CHKERRQ(ierr); 181 ierr = VecDuplicate(x,&s1);CHKERRQ(ierr); 182 ierr = VecDuplicate(x,&s2);CHKERRQ(ierr); 183 184 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr); 185 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 186 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 187 ierr = VecSet(u,one);CHKERRQ(ierr); 188 189 /* Test MatNorm() */ 190 ierr = MatNorm(A,NORM_FROBENIUS,&r1);CHKERRQ(ierr); 191 ierr = MatNorm(sA,NORM_FROBENIUS,&r2);CHKERRQ(ierr); 192 rnorm = PetscAbsReal(r1-r2)/r2; 193 if (rnorm > tol && !rank) { 194 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr); 195 } 196 ierr = MatNorm(A,NORM_INFINITY,&r1);CHKERRQ(ierr); 197 ierr = MatNorm(sA,NORM_INFINITY,&r2);CHKERRQ(ierr); 198 rnorm = PetscAbsReal(r1-r2)/r2; 199 if (rnorm > tol && !rank) { 200 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr); 201 } 202 ierr = MatNorm(A,NORM_1,&r1);CHKERRQ(ierr); 203 ierr = MatNorm(sA,NORM_1,&r2);CHKERRQ(ierr); 204 rnorm = PetscAbsReal(r1-r2)/r2; 205 if (rnorm > tol && !rank) { 206 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);CHKERRQ(ierr); 207 } 208 209 /* Test MatGetOwnershipRange() */ 210 ierr = MatGetOwnershipRange(sA,&rstart,&rend);CHKERRQ(ierr); 211 ierr = MatGetOwnershipRange(A,&i2,&j2);CHKERRQ(ierr); 212 i2 -= rstart; j2 -= rend; 213 if (i2 || j2) { 214 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);CHKERRQ(ierr); 215 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 216 } 217 218 /* Test MatDiagonalScale() */ 219 ierr = MatDiagonalScale(A,x,x);CHKERRQ(ierr); 220 ierr = MatDiagonalScale(sA,x,x);CHKERRQ(ierr); 221 ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr); 222 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale"); 223 224 /* Test MatGetDiagonal(), MatScale() */ 225 ierr = MatGetDiagonal(A,s1);CHKERRQ(ierr); 226 ierr = MatGetDiagonal(sA,s2);CHKERRQ(ierr); 227 ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr); 228 ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr); 229 r1 -= r2; 230 if (r1<-tol || r1>tol) { 231 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);CHKERRQ(ierr); 232 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 233 } 234 235 ierr = MatScale(A,alpha);CHKERRQ(ierr); 236 ierr = MatScale(sA,alpha);CHKERRQ(ierr); 237 238 /* Test MatGetRowMaxAbs() */ 239 ierr = MatGetRowMaxAbs(A,s1,NULL);CHKERRQ(ierr); 240 ierr = MatGetRowMaxAbs(sA,s2,NULL);CHKERRQ(ierr); 241 242 ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr); 243 ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr); 244 r1 -= r2; 245 if (r1<-tol || r1>tol) { 246 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");CHKERRQ(ierr); 247 } 248 249 /* Test MatMult(), MatMultAdd() */ 250 ierr = MatMultEqual(A,sA,10,&flg);CHKERRQ(ierr); 251 if (!flg) { 252 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);CHKERRQ(ierr); 253 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 254 } 255 256 ierr = MatMultAddEqual(A,sA,10,&flg);CHKERRQ(ierr); 257 if (!flg) { 258 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);CHKERRQ(ierr); 259 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 260 } 261 262 /* Test MatMultTranspose(), MatMultTransposeAdd() */ 263 for (i=0; i<10; i++) { 264 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 265 ierr = MatMultTranspose(A,x,s1);CHKERRQ(ierr); 266 ierr = MatMultTranspose(sA,x,s2);CHKERRQ(ierr); 267 ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr); 268 ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr); 269 r1 -= r2; 270 if (r1<-tol || r1>tol) { 271 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);CHKERRQ(ierr); 272 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 273 } 274 } 275 for (i=0; i<10; i++) { 276 ierr = VecSetRandom(x,rctx);CHKERRQ(ierr); 277 ierr = VecSetRandom(y,rctx);CHKERRQ(ierr); 278 ierr = MatMultTransposeAdd(A,x,y,s1);CHKERRQ(ierr); 279 ierr = MatMultTransposeAdd(sA,x,y,s2);CHKERRQ(ierr); 280 ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr); 281 ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr); 282 r1 -= r2; 283 if (r1<-tol || r1>tol) { 284 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);CHKERRQ(ierr); 285 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 286 } 287 } 288 289 /* Test MatDuplicate() */ 290 ierr = MatDuplicate(sA,MAT_COPY_VALUES,&sB);CHKERRQ(ierr); 291 ierr = MatEqual(sA,sB,&flg);CHKERRQ(ierr); 292 if (!flg) { 293 ierr = PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");CHKERRQ(ierr); 294 } 295 ierr = MatMultEqual(sA,sB,5,&flg);CHKERRQ(ierr); 296 if (!flg) { 297 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);CHKERRQ(ierr); 298 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 299 } 300 ierr = MatMultAddEqual(sA,sB,5,&flg);CHKERRQ(ierr); 301 if (!flg) { 302 ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);CHKERRQ(ierr); 303 ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr); 304 } 305 ierr = MatDestroy(&sB);CHKERRQ(ierr); 306 ierr = VecDestroy(&u);CHKERRQ(ierr); 307 ierr = VecDestroy(&x);CHKERRQ(ierr); 308 ierr = VecDestroy(&y);CHKERRQ(ierr); 309 ierr = VecDestroy(&s1);CHKERRQ(ierr); 310 ierr = VecDestroy(&s2);CHKERRQ(ierr); 311 ierr = MatDestroy(&sA);CHKERRQ(ierr); 312 ierr = MatDestroy(&A);CHKERRQ(ierr); 313 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 314 315 ierr = PetscFinalize(); 316 return ierr; 317 } 318 319 /*TEST 320 321 test: 322 nsize: {{1 3}} 323 args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular 324 325 TEST*/ 326