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