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=1; 14 PetscMPIInt size,rank; 15 PetscBool flg; 16 17 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 18 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL)); 19 PetscCall(PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL)); 20 PetscCall(PetscOptionsGetInt(NULL,NULL,"-prob",&prob,NULL)); 21 22 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 24 25 /* Create a BAIJ matrix A */ 26 n = mbs*bs; 27 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 28 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n)); 29 PetscCall(MatSetType(A,MATBAIJ)); 30 PetscCall(MatSetFromOptions(A)); 31 PetscCall(MatMPIBAIJSetPreallocation(A,bs,d_nz,NULL,o_nz,NULL)); 32 PetscCall(MatSeqBAIJSetPreallocation(A,bs,d_nz,NULL)); 33 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 34 35 if (bs == 1) { 36 if (prob == 1) { /* tridiagonal matrix */ 37 value[0] = -1.0; value[1] = 0.0; value[2] = -1.0; 38 for (i=1; i<n-1; i++) { 39 col[0] = i-1; col[1] = i; col[2] = i+1; 40 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 41 } 42 i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1; 43 value[0]= 0.1; value[1]=-1.0; value[2]=0.0; 44 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 45 46 i = 0; col[0] = 0; col[1] = 1; col[2]=n-1; 47 value[0] = 0.0; value[1] = -1.0; value[2]=0.1; 48 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 49 } else if (prob ==2) { /* matrix for the five point stencil */ 50 n1 = (int) PetscSqrtReal((PetscReal)n); 51 for (i=0; i<n1; i++) { 52 for (j=0; j<n1; j++) { 53 Ii = j + n1*i; 54 if (i>0) {J = Ii - n1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES));} 55 if (i<n1-1) {J = Ii + n1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES));} 56 if (j>0) {J = Ii - 1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES));} 57 if (j<n1-1) {J = Ii + 1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES));} 58 PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES)); 59 } 60 } 61 } 62 /* end of if (bs == 1) */ 63 } else { /* bs > 1 */ 64 for (block=0; block<n/bs; block++) { 65 /* diagonal blocks */ 66 value[0] = -1.0; value[1] = 4.0; value[2] = -1.0; 67 for (i=1+block*bs; i<bs-1+block*bs; i++) { 68 col[0] = i-1; col[1] = i; col[2] = i+1; 69 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 70 } 71 i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs; 72 value[0]=-1.0; value[1]=4.0; 73 PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 74 75 i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs; 76 value[0]=4.0; value[1] = -1.0; 77 PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 78 } 79 /* off-diagonal blocks */ 80 value[0]=-1.0; 81 for (i=0; i<(n/bs-1)*bs; i++) { 82 col[0]=i+bs; 83 PetscCall(MatSetValues(A,1,&i,1,col,value,INSERT_VALUES)); 84 col[0]=i; row=i+bs; 85 PetscCall(MatSetValues(A,1,&row,1,col,value,INSERT_VALUES)); 86 } 87 } 88 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 89 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 90 PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 91 92 /* Get SBAIJ matrix sA from A */ 93 PetscCall(MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA)); 94 95 /* Test MatGetSize(), MatGetLocalSize() */ 96 PetscCall(MatGetSize(sA, &i,&j)); 97 PetscCall(MatGetSize(A, &i2,&j2)); 98 i -= i2; j -= j2; 99 if (i || j) { 100 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank)); 101 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 102 } 103 104 PetscCall(MatGetLocalSize(sA, &i,&j)); 105 PetscCall(MatGetLocalSize(A, &i2,&j2)); 106 i2 -= i; j2 -= j; 107 if (i2 || j2) { 108 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank)); 109 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 110 } 111 112 /* vectors */ 113 /*--------------------*/ 114 /* i is obtained from MatGetLocalSize() */ 115 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 116 PetscCall(VecSetSizes(x,i,PETSC_DECIDE)); 117 PetscCall(VecSetFromOptions(x)); 118 PetscCall(VecDuplicate(x,&y)); 119 PetscCall(VecDuplicate(x,&u)); 120 PetscCall(VecDuplicate(x,&s1)); 121 PetscCall(VecDuplicate(x,&s2)); 122 123 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rctx)); 124 PetscCall(PetscRandomSetFromOptions(rctx)); 125 PetscCall(VecSetRandom(x,rctx)); 126 PetscCall(VecSet(u,one)); 127 128 /* Test MatNorm() */ 129 PetscCall(MatNorm(A,NORM_FROBENIUS,&r1)); 130 PetscCall(MatNorm(sA,NORM_FROBENIUS,&r2)); 131 rnorm = PetscAbsReal(r1-r2)/r2; 132 if (rnorm > tol && rank == 0) { 133 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs)); 134 } 135 PetscCall(MatNorm(A,NORM_INFINITY,&r1)); 136 PetscCall(MatNorm(sA,NORM_INFINITY,&r2)); 137 rnorm = PetscAbsReal(r1-r2)/r2; 138 if (rnorm > tol && rank == 0) { 139 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs)); 140 } 141 PetscCall(MatNorm(A,NORM_1,&r1)); 142 PetscCall(MatNorm(sA,NORM_1,&r2)); 143 rnorm = PetscAbsReal(r1-r2)/r2; 144 if (rnorm > tol && rank == 0) { 145 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n",(double)r1,(double)r2,bs)); 146 } 147 148 /* Test MatGetOwnershipRange() */ 149 PetscCall(MatGetOwnershipRange(sA,&rstart,&rend)); 150 PetscCall(MatGetOwnershipRange(A,&i2,&j2)); 151 i2 -= rstart; j2 -= rend; 152 if (i2 || j2) { 153 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank)); 154 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 155 } 156 157 /* Test MatDiagonalScale() */ 158 PetscCall(MatDiagonalScale(A,x,x)); 159 PetscCall(MatDiagonalScale(sA,x,x)); 160 PetscCall(MatMultEqual(A,sA,10,&flg)); 161 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale"); 162 163 /* Test MatGetDiagonal(), MatScale() */ 164 PetscCall(MatGetDiagonal(A,s1)); 165 PetscCall(MatGetDiagonal(sA,s2)); 166 PetscCall(VecNorm(s1,NORM_1,&r1)); 167 PetscCall(VecNorm(s2,NORM_1,&r2)); 168 r1 -= r2; 169 if (r1<-tol || r1>tol) { 170 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1)); 171 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 172 } 173 174 PetscCall(MatScale(A,alpha)); 175 PetscCall(MatScale(sA,alpha)); 176 177 /* Test MatGetRowMaxAbs() */ 178 PetscCall(MatGetRowMaxAbs(A,s1,NULL)); 179 PetscCall(MatGetRowMaxAbs(sA,s2,NULL)); 180 181 PetscCall(VecNorm(s1,NORM_1,&r1)); 182 PetscCall(VecNorm(s2,NORM_1,&r2)); 183 r1 -= r2; 184 if (r1<-tol || r1>tol) { 185 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n")); 186 } 187 188 /* Test MatMult(), MatMultAdd() */ 189 PetscCall(MatMultEqual(A,sA,10,&flg)); 190 if (!flg) { 191 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank)); 192 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 193 } 194 195 PetscCall(MatMultAddEqual(A,sA,10,&flg)); 196 if (!flg) { 197 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank)); 198 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 199 } 200 201 /* Test MatMultTranspose(), MatMultTransposeAdd() */ 202 for (i=0; i<10; i++) { 203 PetscCall(VecSetRandom(x,rctx)); 204 PetscCall(MatMultTranspose(A,x,s1)); 205 PetscCall(MatMultTranspose(sA,x,s2)); 206 PetscCall(VecNorm(s1,NORM_1,&r1)); 207 PetscCall(VecNorm(s2,NORM_1,&r2)); 208 r1 -= r2; 209 if (r1<-tol || r1>tol) { 210 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1)); 211 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 212 } 213 } 214 for (i=0; i<10; i++) { 215 PetscCall(VecSetRandom(x,rctx)); 216 PetscCall(VecSetRandom(y,rctx)); 217 PetscCall(MatMultTransposeAdd(A,x,y,s1)); 218 PetscCall(MatMultTransposeAdd(sA,x,y,s2)); 219 PetscCall(VecNorm(s1,NORM_1,&r1)); 220 PetscCall(VecNorm(s2,NORM_1,&r2)); 221 r1 -= r2; 222 if (r1<-tol || r1>tol) { 223 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1)); 224 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 225 } 226 } 227 228 /* Test MatDuplicate() */ 229 PetscCall(MatDuplicate(sA,MAT_COPY_VALUES,&sB)); 230 PetscCall(MatEqual(sA,sB,&flg)); 231 if (!flg) { 232 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n")); 233 } 234 PetscCall(MatMultEqual(sA,sB,5,&flg)); 235 if (!flg) { 236 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank)); 237 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 238 } 239 PetscCall(MatMultAddEqual(sA,sB,5,&flg)); 240 if (!flg) { 241 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank)); 242 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 243 } 244 PetscCall(MatDestroy(&sB)); 245 PetscCall(VecDestroy(&u)); 246 PetscCall(VecDestroy(&x)); 247 PetscCall(VecDestroy(&y)); 248 PetscCall(VecDestroy(&s1)); 249 PetscCall(VecDestroy(&s2)); 250 PetscCall(MatDestroy(&sA)); 251 PetscCall(MatDestroy(&A)); 252 PetscCall(PetscRandomDestroy(&rctx)); 253 PetscCall(PetscFinalize()); 254 return 0; 255 } 256 257 /*TEST 258 259 test: 260 nsize: {{1 3}} 261 args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular -prob {{1 2}} 262 263 TEST*/ 264