1 #include <petscmat.h> 2 3 static char help[] = "Tests MatMat operations with MAT_REUSE_MATRIX and already allocated dense result.\n\n"; 4 5 static PetscScalar MAGIC_NUMBER = 12345; 6 7 static PetscErrorCode CheckLocal(Mat A, Mat B, PetscScalar *a, PetscScalar *b) 8 { 9 PetscBool wA = PETSC_FALSE, wB = PETSC_FALSE; 10 PetscBool wAv = PETSC_FALSE, wBv = PETSC_FALSE; 11 PetscInt lda,i,j,m,n; 12 13 PetscFunctionBegin; 14 if (a) { 15 const PetscScalar *Aa; 16 CHKERRQ(MatDenseGetArrayRead(A,&Aa)); 17 wA = (PetscBool)(a != Aa); 18 CHKERRQ(MatDenseGetLDA(A,&lda)); 19 CHKERRQ(MatGetLocalSize(A,&m,&n)); 20 for (j=0;j<n;j++) { 21 for (i=m;i<lda;i++) { 22 if (Aa[j*lda +i] != MAGIC_NUMBER) wAv = PETSC_TRUE; 23 } 24 } 25 CHKERRQ(MatDenseRestoreArrayRead(A,&Aa)); 26 } 27 if (b) { 28 const PetscScalar *Bb; 29 CHKERRQ(MatDenseGetArrayRead(B,&Bb)); 30 wB = (PetscBool)(b != Bb); 31 CHKERRQ(MatDenseGetLDA(B,&lda)); 32 CHKERRQ(MatGetLocalSize(B,&m,&n)); 33 for (j=0;j<n;j++) { 34 for (i=m;i<lda;i++) { 35 if (Bb[j*lda +i] != MAGIC_NUMBER) wBv = PETSC_TRUE; 36 } 37 } 38 CHKERRQ(MatDenseRestoreArrayRead(B,&Bb)); 39 } 40 PetscCheckFalse(wA || wB,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong array in first Mat? %d, Wrong array in second Mat? %d",wA,wB); 41 PetscCheckFalse(wAv || wBv,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong data in first Mat? %d, Wrong data in second Mat? %d",wAv,wBv); 42 PetscFunctionReturn(0); 43 } 44 45 typedef struct { 46 Mat A; 47 Mat P; 48 Mat R; 49 } proj_data; 50 51 PetscErrorCode proj_destroy(void *ctx) 52 { 53 proj_data *userdata = (proj_data*)ctx; 54 55 PetscFunctionBegin; 56 PetscCheck(userdata,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing userdata"); 57 CHKERRQ(MatDestroy(&userdata->A)); 58 CHKERRQ(MatDestroy(&userdata->P)); 59 CHKERRQ(MatDestroy(&userdata->R)); 60 CHKERRQ(PetscFree(userdata)); 61 PetscFunctionReturn(0); 62 } 63 64 PetscErrorCode proj_mult(Mat S, Vec X, Vec Y) 65 { 66 Mat A,R,P; 67 Vec Ax,Ay; 68 Vec Px,Py; 69 proj_data *userdata; 70 71 PetscFunctionBegin; 72 CHKERRQ(MatShellGetContext(S,&userdata)); 73 PetscCheck(userdata,PetscObjectComm((PetscObject)S),PETSC_ERR_PLIB,"Missing userdata"); 74 A = userdata->A; 75 R = userdata->R; 76 P = userdata->P; 77 PetscCheck(A,PetscObjectComm((PetscObject)S),PETSC_ERR_PLIB,"Missing matrix"); 78 PetscCheckFalse(!R && !P,PetscObjectComm((PetscObject)S),PETSC_ERR_PLIB,"Missing projectors"); 79 PetscCheckFalse(R && P,PetscObjectComm((PetscObject)S),PETSC_ERR_PLIB,"Both projectors"); 80 CHKERRQ(MatCreateVecs(A,&Ax,&Ay)); 81 if (R) { 82 CHKERRQ(MatCreateVecs(R,&Py,&Px)); 83 } else { 84 CHKERRQ(MatCreateVecs(P,&Px,&Py)); 85 } 86 CHKERRQ(VecCopy(X,Px)); 87 if (P) { 88 CHKERRQ(MatMult(P,Px,Py)); 89 } else { 90 CHKERRQ(MatMultTranspose(R,Px,Py)); 91 } 92 CHKERRQ(VecCopy(Py,Ax)); 93 CHKERRQ(MatMult(A,Ax,Ay)); 94 CHKERRQ(VecCopy(Ay,Py)); 95 if (P) { 96 CHKERRQ(MatMultTranspose(P,Py,Px)); 97 } else { 98 CHKERRQ(MatMult(R,Py,Px)); 99 } 100 CHKERRQ(VecCopy(Px,Y)); 101 CHKERRQ(VecDestroy(&Px)); 102 CHKERRQ(VecDestroy(&Py)); 103 CHKERRQ(VecDestroy(&Ax)); 104 CHKERRQ(VecDestroy(&Ay)); 105 PetscFunctionReturn(0); 106 } 107 108 PetscErrorCode MyPtShellPMultSymbolic(Mat S, Mat P, Mat PtAP, void** ctx) 109 { 110 proj_data *userdata; 111 112 PetscFunctionBegin; 113 CHKERRQ(PetscNew(&userdata)); 114 CHKERRQ(MatShellSetContext(PtAP,userdata)); 115 *ctx = (void *)userdata; 116 PetscFunctionReturn(0); 117 } 118 119 PetscErrorCode MyPtShellPMultNumeric(Mat S, Mat P, Mat PtAP, void *ctx) 120 { 121 Mat A; 122 proj_data *userdata = (proj_data*)ctx; 123 124 PetscFunctionBegin; 125 CHKERRQ(MatShellGetContext(S,&A)); 126 CHKERRQ(PetscObjectReference((PetscObject)A)); 127 CHKERRQ(PetscObjectReference((PetscObject)P)); 128 CHKERRQ(MatDestroy(&userdata->A)); 129 CHKERRQ(MatDestroy(&userdata->P)); 130 CHKERRQ(MatDestroy(&userdata->R)); 131 userdata->A = A; 132 userdata->P = P; 133 CHKERRQ(MatShellSetOperation(PtAP,MATOP_MULT,(void (*)(void))proj_mult)); 134 CHKERRQ(MatSetUp(PtAP)); 135 CHKERRQ(MatAssemblyBegin(PtAP,MAT_FINAL_ASSEMBLY)); 136 CHKERRQ(MatAssemblyEnd(PtAP,MAT_FINAL_ASSEMBLY)); 137 PetscFunctionReturn(0); 138 } 139 140 PetscErrorCode MyRShellRtMultSymbolic(Mat S, Mat R, Mat RARt, void **ctx) 141 { 142 proj_data *userdata; 143 144 PetscFunctionBegin; 145 CHKERRQ(PetscNew(&userdata)); 146 CHKERRQ(MatShellSetContext(RARt,userdata)); 147 *ctx = (void *)userdata; 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode MyRShellRtMultNumeric(Mat S, Mat R, Mat RARt, void *ctx) 152 { 153 Mat A; 154 proj_data *userdata = (proj_data*)ctx; 155 156 PetscFunctionBegin; 157 CHKERRQ(MatShellGetContext(S,&A)); 158 CHKERRQ(PetscObjectReference((PetscObject)A)); 159 CHKERRQ(PetscObjectReference((PetscObject)R)); 160 CHKERRQ(MatDestroy(&userdata->A)); 161 CHKERRQ(MatDestroy(&userdata->P)); 162 CHKERRQ(MatDestroy(&userdata->R)); 163 userdata->A = A; 164 userdata->R = R; 165 CHKERRQ(MatShellSetOperation(RARt,MATOP_MULT,(void (*)(void))proj_mult)); 166 CHKERRQ(MatSetUp(RARt)); 167 CHKERRQ(MatAssemblyBegin(RARt,MAT_FINAL_ASSEMBLY)); 168 CHKERRQ(MatAssemblyEnd(RARt,MAT_FINAL_ASSEMBLY)); 169 PetscFunctionReturn(0); 170 } 171 172 PetscErrorCode MyMatShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx) 173 { 174 Mat A; 175 176 PetscFunctionBegin; 177 CHKERRQ(MatShellGetContext(S,&A)); 178 CHKERRQ(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C)); 179 PetscFunctionReturn(0); 180 } 181 182 PetscErrorCode MyMatTransposeShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx) 183 { 184 Mat A; 185 186 PetscFunctionBegin; 187 CHKERRQ(MatShellGetContext(S,&A)); 188 CHKERRQ(MatTransposeMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C)); 189 PetscFunctionReturn(0); 190 } 191 192 PetscErrorCode MyMatShellMatTransposeMultNumeric(Mat S, Mat B, Mat C, void *ctx) 193 { 194 Mat A; 195 196 PetscFunctionBegin; 197 CHKERRQ(MatShellGetContext(S,&A)); 198 CHKERRQ(MatMatTransposeMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C)); 199 PetscFunctionReturn(0); 200 } 201 202 int main(int argc,char **args) 203 { 204 Mat X,B,A,Bt,T,T2,PtAP = NULL,RARt = NULL, R = NULL; 205 Vec r,l,rs,ls; 206 PetscInt m,n,k,M = 10,N = 10,K = 5, ldx = 3, ldb = 5, ldr = 4; 207 const char *deft = MATAIJ; 208 char mattype[256]; 209 PetscBool flg,symm = PETSC_FALSE,testtt = PETSC_TRUE, testnest = PETSC_TRUE, testtranspose = PETSC_TRUE, testcircular = PETSC_FALSE, local = PETSC_TRUE; 210 PetscBool testhtranspose = PETSC_TRUE; 211 PetscBool xgpu = PETSC_FALSE, bgpu = PETSC_FALSE, testshellops = PETSC_FALSE, testproj = PETSC_TRUE, testrart = PETSC_TRUE, testmatmatt = PETSC_TRUE, testmattmat = PETSC_TRUE; 212 PetscScalar *dataX = NULL,*dataB = NULL, *dataR = NULL, *dataBt = NULL; 213 PetscScalar *aX,*aB,*aBt; 214 PetscReal err; 215 PetscErrorCode ierr; 216 217 ierr = PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr; 218 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 219 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 220 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL)); 221 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL)); 222 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-local",&local,NULL)); 223 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-ldx",&ldx,NULL)); 224 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-ldb",&ldb,NULL)); 225 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-ldr",&ldr,NULL)); 226 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testtranspose",&testtranspose,NULL)); 227 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testnest",&testnest,NULL)); 228 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testtt",&testtt,NULL)); 229 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testcircular",&testcircular,NULL)); 230 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testshellops",&testshellops,NULL)); 231 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testproj",&testproj,NULL)); 232 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testrart",&testrart,NULL)); 233 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testmatmatt",&testmatmatt,NULL)); 234 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-testmattmat",&testmattmat,NULL)); 235 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-xgpu",&xgpu,NULL)); 236 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-bgpu",&bgpu,NULL)); 237 CHKERRQ(PetscOptionsGetScalar(NULL,NULL,"-magic_number",&MAGIC_NUMBER,NULL)); 238 if (M != N) testproj = PETSC_FALSE; 239 240 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 241 CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N)); 242 CHKERRQ(MatSetType(A,MATAIJ)); 243 CHKERRQ(MatSetUp(A)); 244 CHKERRQ(MatSetRandom(A,NULL)); 245 if (M==N && symm) { 246 Mat AT; 247 248 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&AT)); 249 CHKERRQ(MatAXPY(A,1.0,AT,DIFFERENT_NONZERO_PATTERN)); 250 CHKERRQ(MatDestroy(&AT)); 251 CHKERRQ(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 252 } 253 CHKERRQ(MatViewFromOptions(A,NULL,"-A_init_view")); 254 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","","");CHKERRQ(ierr); 255 CHKERRQ(PetscOptionsFList("-A_mat_type","Matrix type","MatSetType",MatList,deft,mattype,256,&flg)); 256 ierr = PetscOptionsEnd();CHKERRQ(ierr); 257 if (flg) { 258 Mat A2; 259 260 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 261 CHKERRQ(MatConvert(A,mattype,MAT_INPLACE_MATRIX,&A)); 262 CHKERRQ(MatMultEqual(A,A2,10,&flg)); 263 if (!flg) { 264 Mat AE,A2E; 265 266 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with convert\n")); 267 CHKERRQ(MatComputeOperator(A,MATDENSE,&AE)); 268 CHKERRQ(MatComputeOperator(A2,MATDENSE,&A2E)); 269 CHKERRQ(MatView(AE,NULL)); 270 CHKERRQ(MatView(A2E,NULL)); 271 CHKERRQ(MatAXPY(A2E,-1.0,A,SAME_NONZERO_PATTERN)); 272 CHKERRQ(MatView(A2E,NULL)); 273 CHKERRQ(MatDestroy(&A2E)); 274 CHKERRQ(MatDestroy(&AE)); 275 } 276 CHKERRQ(MatDestroy(&A2)); 277 } 278 CHKERRQ(MatViewFromOptions(A,NULL,"-A_view")); 279 280 CHKERRQ(MatGetLocalSize(A,&m,&n)); 281 if (local) { 282 PetscInt i; 283 284 CHKERRQ(PetscMalloc1((m+ldx)*K,&dataX)); 285 CHKERRQ(PetscMalloc1((n+ldb)*K,&dataB)); 286 for (i=0;i<(m+ldx)*K;i++) dataX[i] = MAGIC_NUMBER; 287 for (i=0;i<(n+ldb)*K;i++) dataB[i] = MAGIC_NUMBER; 288 } 289 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,K,dataB,&B)); 290 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,dataX,&X)); 291 if (local) { 292 CHKERRQ(MatDenseSetLDA(X,m+ldx)); 293 CHKERRQ(MatDenseSetLDA(B,n+ldb)); 294 } 295 CHKERRQ(MatGetLocalSize(B,NULL,&k)); 296 if (local) { 297 PetscInt i; 298 299 CHKERRQ(PetscMalloc1((k+ldr)*N,&dataBt)); 300 for (i=0;i<(k+ldr)*N;i++) dataBt[i] = MAGIC_NUMBER; 301 } 302 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,k,n,K,N,dataBt,&Bt)); 303 if (local) { 304 CHKERRQ(MatDenseSetLDA(Bt,k+ldr)); 305 } 306 307 /* store pointer to dense data for testing */ 308 CHKERRQ(MatDenseGetArrayRead(B,(const PetscScalar**)&dataB)); 309 CHKERRQ(MatDenseGetArrayRead(X,(const PetscScalar**)&dataX)); 310 CHKERRQ(MatDenseGetArrayRead(Bt,(const PetscScalar**)&dataBt)); 311 aX = dataX; 312 aB = dataB; 313 aBt = dataBt; 314 CHKERRQ(MatDenseRestoreArrayRead(Bt,(const PetscScalar**)&dataBt)); 315 CHKERRQ(MatDenseRestoreArrayRead(B,(const PetscScalar**)&dataB)); 316 CHKERRQ(MatDenseRestoreArrayRead(X,(const PetscScalar**)&dataX)); 317 if (local) { 318 dataX = aX; 319 dataB = aB; 320 dataBt = aBt; 321 } 322 323 CHKERRQ(MatSetRandom(X,NULL)); 324 CHKERRQ(MatSetRandom(B,NULL)); 325 CHKERRQ(MatSetRandom(Bt,NULL)); 326 CHKERRQ(CheckLocal(X,NULL,aX,NULL)); 327 CHKERRQ(CheckLocal(Bt,B,aBt,aB)); 328 329 /* convert to CUDA if needed */ 330 if (bgpu) { 331 CHKERRQ(MatConvert(B,MATDENSECUDA,MAT_INPLACE_MATRIX,&B)); 332 CHKERRQ(MatConvert(Bt,MATDENSECUDA,MAT_INPLACE_MATRIX,&Bt)); 333 } 334 if (xgpu) { 335 CHKERRQ(MatConvert(X,MATDENSECUDA,MAT_INPLACE_MATRIX,&X)); 336 } 337 CHKERRQ(CheckLocal(B,X,aB,aX)); 338 339 /* Test MatDenseGetSubMatrix */ 340 { 341 Mat B2,T3,T4; 342 343 CHKERRQ(MatDuplicate(B,MAT_COPY_VALUES,&B2)); 344 CHKERRQ(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&T4)); 345 CHKERRQ(MatSetRandom(T4,NULL)); 346 CHKERRQ(MatAXPY(B2,1.0,T4,SAME_NONZERO_PATTERN)); 347 CHKERRQ(MatDenseGetSubMatrix(B,PetscMin(1,K),PetscMin(2,K),&T)); 348 CHKERRQ(MatDenseGetSubMatrix(T4,PetscMin(1,K),PetscMin(2,K),&T2)); 349 CHKERRQ(MatDenseGetSubMatrix(B2,PetscMin(1,K),PetscMin(2,K),&T3)); 350 CHKERRQ(MatAXPY(T,1.0,T2,SAME_NONZERO_PATTERN)); 351 CHKERRQ(MatAXPY(T3,-1.0,T,SAME_NONZERO_PATTERN)); 352 CHKERRQ(MatNorm(T3,NORM_FROBENIUS,&err)); 353 if (err > PETSC_SMALL) { 354 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with MatDenseGetSubMatrix\n")); 355 CHKERRQ(MatView(T3,NULL)); 356 } 357 CHKERRQ(MatDenseRestoreSubMatrix(B,&T)); 358 CHKERRQ(MatDenseRestoreSubMatrix(T4,&T2)); 359 CHKERRQ(MatDenseRestoreSubMatrix(B2,&T3)); 360 CHKERRQ(CheckLocal(B,NULL,aB,NULL)); 361 CHKERRQ(MatDestroy(&B2)); 362 CHKERRQ(MatDestroy(&T4)); 363 } 364 365 /* Test reusing a previously allocated dense buffer */ 366 CHKERRQ(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X)); 367 CHKERRQ(CheckLocal(B,X,aB,aX)); 368 CHKERRQ(MatMatMultEqual(A,B,X,10,&flg)); 369 if (!flg) { 370 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage\n")); 371 CHKERRQ(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 372 CHKERRQ(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN)); 373 CHKERRQ(MatView(T,NULL)); 374 CHKERRQ(MatDestroy(&T)); 375 } 376 377 /* Test MatTransposeMat and MatMatTranspose */ 378 if (testmattmat) { 379 CHKERRQ(MatTransposeMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B)); 380 CHKERRQ(CheckLocal(B,X,aB,aX)); 381 CHKERRQ(MatTransposeMatMultEqual(A,X,B,10,&flg)); 382 if (!flg) { 383 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatTransposeMat)\n")); 384 CHKERRQ(MatTransposeMatMult(A,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B)); 385 CHKERRQ(MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN)); 386 CHKERRQ(MatView(T,NULL)); 387 CHKERRQ(MatDestroy(&T)); 388 } 389 } 390 if (testmatmatt) { 391 CHKERRQ(MatMatTransposeMult(A,Bt,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X)); 392 CHKERRQ(CheckLocal(Bt,X,aBt,aX)); 393 CHKERRQ(MatMatTransposeMultEqual(A,Bt,X,10,&flg)); 394 if (!flg) { 395 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatMatTranspose)\n")); 396 CHKERRQ(MatMatTransposeMult(A,Bt,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 397 CHKERRQ(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN)); 398 CHKERRQ(MatView(T,NULL)); 399 CHKERRQ(MatDestroy(&T)); 400 } 401 } 402 403 /* Test projection operations (PtAP and RARt) */ 404 if (testproj) { 405 CHKERRQ(MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&PtAP)); 406 CHKERRQ(MatPtAPMultEqual(A,B,PtAP,10,&flg)); 407 if (!flg) { 408 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with PtAP\n")); 409 CHKERRQ(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 410 CHKERRQ(MatTransposeMatMult(B,T,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T2)); 411 CHKERRQ(MatAXPY(T2,-1.0,PtAP,SAME_NONZERO_PATTERN)); 412 CHKERRQ(MatView(T2,NULL)); 413 CHKERRQ(MatDestroy(&T2)); 414 CHKERRQ(MatDestroy(&T)); 415 } 416 CHKERRQ(PetscMalloc1((k+ldr)*M,&dataR)); 417 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,m,K,M,dataR,&R)); 418 CHKERRQ(MatDenseSetLDA(R,k+ldr)); 419 CHKERRQ(MatSetRandom(R,NULL)); 420 if (testrart) { /* fails for AIJCUSPARSE because RA operation is not defined */ 421 CHKERRQ(MatRARt(A,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&RARt)); 422 CHKERRQ(MatRARtMultEqual(A,R,RARt,10,&flg)); 423 if (!flg) { 424 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with RARt\n")); 425 CHKERRQ(MatMatTransposeMult(A,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 426 CHKERRQ(MatMatMult(R,T,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T2)); 427 CHKERRQ(MatAXPY(T2,-1.0,RARt,SAME_NONZERO_PATTERN)); 428 CHKERRQ(MatView(T2,NULL)); 429 CHKERRQ(MatDestroy(&T2)); 430 CHKERRQ(MatDestroy(&T)); 431 } 432 } 433 } 434 435 /* Test MatDenseGetColumnVec and friends */ 436 CHKERRQ(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X)); 437 CHKERRQ(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 438 CHKERRQ(MatDuplicate(T,MAT_DO_NOT_COPY_VALUES,&T2)); 439 for (k=0;k<K;k++) { 440 Vec Xv,Tv,T2v; 441 442 CHKERRQ(MatDenseGetColumnVecRead(X,k,&Xv)); 443 CHKERRQ(MatDenseGetColumnVec(T,k,&Tv)); 444 CHKERRQ(MatDenseGetColumnVecWrite(T2,k,&T2v)); 445 CHKERRQ(VecCopy(Xv,T2v)); 446 CHKERRQ(VecAXPY(Tv,-1.,Xv)); 447 CHKERRQ(MatDenseRestoreColumnVecRead(X,k,&Xv)); 448 CHKERRQ(MatDenseRestoreColumnVec(T,k,&Tv)); 449 CHKERRQ(MatDenseRestoreColumnVecWrite(T2,k,&T2v)); 450 } 451 CHKERRQ(MatNorm(T,NORM_FROBENIUS,&err)); 452 if (err > PETSC_SMALL) { 453 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with MatDenseGetColumnVec\n")); 454 CHKERRQ(MatView(T,NULL)); 455 } 456 CHKERRQ(MatAXPY(T2,-1.,X,SAME_NONZERO_PATTERN)); 457 CHKERRQ(MatNorm(T2,NORM_FROBENIUS,&err)); 458 if (err > PETSC_SMALL) { 459 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with MatDenseGetColumnVecWrite\n")); 460 CHKERRQ(MatView(T2,NULL)); 461 } 462 CHKERRQ(MatDestroy(&T)); 463 CHKERRQ(MatDestroy(&T2)); 464 465 /* Test with MatShell */ 466 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&T)); 467 CHKERRQ(MatConvert(T,MATSHELL,MAT_INITIAL_MATRIX,&T2)); 468 CHKERRQ(MatDestroy(&T)); 469 470 /* scale matrix */ 471 CHKERRQ(MatScale(A,2.0)); 472 CHKERRQ(MatScale(T2,2.0)); 473 CHKERRQ(MatCreateVecs(A,&r,&l)); 474 CHKERRQ(VecSetRandom(r,NULL)); 475 CHKERRQ(VecSetRandom(l,NULL)); 476 CHKERRQ(MatCreateVecs(T2,&rs,&ls)); 477 CHKERRQ(VecCopy(r,rs)); 478 CHKERRQ(VecCopy(l,ls)); 479 if (testproj) { 480 CHKERRQ(MatDiagonalScale(A,r,r)); 481 CHKERRQ(MatDiagonalScale(T2,rs,rs)); 482 } else { 483 CHKERRQ(MatDiagonalScale(A,l,r)); 484 CHKERRQ(MatDiagonalScale(T2,ls,rs)); 485 } 486 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&T)); 487 CHKERRQ(MatAXPY(A,4.5,T,SAME_NONZERO_PATTERN)); 488 CHKERRQ(MatAXPY(T2,4.5,T,DIFFERENT_NONZERO_PATTERN)); 489 CHKERRQ(MatMultEqual(T2,A,10,&flg)); 490 if (!flg) { 491 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with MATSHELL (MatMult)\n")); 492 } 493 CHKERRQ(MatMultTransposeEqual(T2,A,10,&flg)); 494 if (!flg) { 495 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with MATSHELL (MatMultTranspose)\n")); 496 } 497 CHKERRQ(MatDestroy(&T)); 498 CHKERRQ(VecDestroy(&ls)); 499 CHKERRQ(VecDestroy(&rs)); 500 CHKERRQ(VecDestroy(&l)); 501 CHKERRQ(VecDestroy(&r)); 502 503 /* recompute projections, test reusage */ 504 if (PtAP) CHKERRQ(MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&PtAP)); 505 if (RARt) CHKERRQ(MatRARt(A,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&RARt)); 506 if (testshellops) { /* test callbacks for user defined MatProducts */ 507 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_AB,NULL,MyMatShellMatMultNumeric,NULL,MATDENSE,MATDENSE)); 508 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_AB,NULL,MyMatShellMatMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA)); 509 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_AtB,NULL,MyMatTransposeShellMatMultNumeric,NULL,MATDENSE,MATDENSE)); 510 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_AtB,NULL,MyMatTransposeShellMatMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA)); 511 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_ABt,NULL,MyMatShellMatTransposeMultNumeric,NULL,MATDENSE,MATDENSE)); 512 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_ABt,NULL,MyMatShellMatTransposeMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA)); 513 if (testproj) { 514 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_PtAP,MyPtShellPMultSymbolic,MyPtShellPMultNumeric,proj_destroy,MATDENSE,MATSHELL)); 515 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_PtAP,MyPtShellPMultSymbolic,MyPtShellPMultNumeric,proj_destroy,MATDENSECUDA,MATSHELL)); 516 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_RARt,MyRShellRtMultSymbolic,MyRShellRtMultNumeric,proj_destroy,MATDENSE,MATSHELL)); 517 CHKERRQ(MatShellSetMatProductOperation(T2,MATPRODUCT_RARt,MyRShellRtMultSymbolic,MyRShellRtMultNumeric,proj_destroy,MATDENSECUDA,MATSHELL)); 518 } 519 } 520 CHKERRQ(CheckLocal(B,X,aB,aX)); 521 /* we either use the shell operations or the loop over columns code, applying the operator */ 522 CHKERRQ(MatMatMult(T2,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X)); 523 CHKERRQ(CheckLocal(B,X,aB,aX)); 524 CHKERRQ(MatMatMultEqual(T2,B,X,10,&flg)); 525 if (!flg) { 526 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MATSHELL)\n")); 527 CHKERRQ(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 528 CHKERRQ(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN)); 529 CHKERRQ(MatView(T,NULL)); 530 CHKERRQ(MatDestroy(&T)); 531 } 532 if (testproj) { 533 CHKERRQ(MatPtAPMultEqual(T2,B,PtAP,10,&flg)); 534 if (!flg) { 535 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with PtAP (MATSHELL)\n")); 536 } 537 if (testshellops) { /* projections fail if the product operations are not specified */ 538 CHKERRQ(MatPtAP(T2,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 539 CHKERRQ(MatPtAP(T2,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&T)); 540 CHKERRQ(MatPtAPMultEqual(T2,B,T,10,&flg)); 541 if (!flg) { 542 Mat TE; 543 544 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with PtAP (MATSHELL user defined)\n")); 545 CHKERRQ(MatComputeOperator(T,MATDENSE,&TE)); 546 CHKERRQ(MatView(TE,NULL)); 547 CHKERRQ(MatView(PtAP,NULL)); 548 CHKERRQ(MatAXPY(TE,-1.0,PtAP,SAME_NONZERO_PATTERN)); 549 CHKERRQ(MatView(TE,NULL)); 550 CHKERRQ(MatDestroy(&TE)); 551 } 552 CHKERRQ(MatDestroy(&T)); 553 } 554 if (RARt) { 555 CHKERRQ(MatRARtMultEqual(T2,R,RARt,10,&flg)); 556 if (!flg) { 557 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with RARt (MATSHELL)\n")); 558 } 559 } 560 if (testshellops) { 561 CHKERRQ(MatRARt(T2,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 562 CHKERRQ(MatRARt(T2,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&T)); 563 CHKERRQ(MatRARtMultEqual(T2,R,T,10,&flg)); 564 if (!flg) { 565 Mat TE; 566 567 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with RARt (MATSHELL user defined)\n")); 568 CHKERRQ(MatComputeOperator(T,MATDENSE,&TE)); 569 CHKERRQ(MatView(TE,NULL)); 570 if (RARt) { 571 CHKERRQ(MatView(RARt,NULL)); 572 CHKERRQ(MatAXPY(TE,-1.0,RARt,SAME_NONZERO_PATTERN)); 573 CHKERRQ(MatView(TE,NULL)); 574 } 575 CHKERRQ(MatDestroy(&TE)); 576 } 577 CHKERRQ(MatDestroy(&T)); 578 } 579 } 580 581 if (testmattmat) { /* we either use the shell operations or the loop over columns code applying the transposed operator */ 582 CHKERRQ(MatTransposeMatMult(T2,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B)); 583 CHKERRQ(CheckLocal(B,X,aB,aX)); 584 CHKERRQ(MatTransposeMatMultEqual(T2,X,B,10,&flg)); 585 if (!flg) { 586 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatTranspose, MATSHELL)\n")); 587 CHKERRQ(MatTransposeMatMult(A,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 588 CHKERRQ(MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN)); 589 CHKERRQ(MatView(T,NULL)); 590 CHKERRQ(MatDestroy(&T)); 591 } 592 } 593 if (testmatmatt && testshellops) { /* only when shell operations are set */ 594 CHKERRQ(MatMatTransposeMult(T2,Bt,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X)); 595 CHKERRQ(CheckLocal(Bt,X,aBt,aX)); 596 CHKERRQ(MatMatTransposeMultEqual(T2,Bt,X,10,&flg)); 597 if (!flg) { 598 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatMatTranspose, MATSHELL)\n")); 599 CHKERRQ(MatMatTransposeMult(A,Bt,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 600 CHKERRQ(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN)); 601 CHKERRQ(MatView(T,NULL)); 602 CHKERRQ(MatDestroy(&T)); 603 } 604 } 605 CHKERRQ(MatDestroy(&T2)); 606 607 if (testnest) { /* test with MatNest */ 608 Mat NA; 609 610 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,1,NULL,1,NULL,&A,&NA)); 611 CHKERRQ(MatViewFromOptions(NA,NULL,"-NA_view")); 612 CHKERRQ(MatMatMult(NA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X)); 613 CHKERRQ(CheckLocal(B,X,aB,aX)); 614 CHKERRQ(MatMatMultEqual(NA,B,X,10,&flg)); 615 if (!flg) { 616 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with Nest\n")); 617 CHKERRQ(MatMatMult(NA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 618 CHKERRQ(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN)); 619 CHKERRQ(MatView(T,NULL)); 620 CHKERRQ(MatDestroy(&T)); 621 } 622 CHKERRQ(MatDestroy(&NA)); 623 } 624 625 if (testtranspose) { /* test with Transpose */ 626 Mat TA; 627 628 CHKERRQ(MatCreateTranspose(A,&TA)); 629 CHKERRQ(MatMatMult(TA,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B)); 630 CHKERRQ(CheckLocal(B,X,aB,aX)); 631 CHKERRQ(MatMatMultEqual(TA,X,B,10,&flg)); 632 if (!flg) { 633 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose\n")); 634 CHKERRQ(MatMatMult(TA,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 635 CHKERRQ(MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN)); 636 CHKERRQ(MatView(T,NULL)); 637 CHKERRQ(MatDestroy(&T)); 638 } 639 CHKERRQ(MatDestroy(&TA)); 640 } 641 642 if (testhtranspose) { /* test with Hermitian Transpose */ 643 Mat TA; 644 645 CHKERRQ(MatCreateHermitianTranspose(A,&TA)); 646 CHKERRQ(MatMatMult(TA,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B)); 647 CHKERRQ(CheckLocal(B,X,aB,aX)); 648 CHKERRQ(MatMatMultEqual(TA,X,B,10,&flg)); 649 if (!flg) { 650 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose\n")); 651 CHKERRQ(MatMatMult(TA,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 652 CHKERRQ(MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN)); 653 CHKERRQ(MatView(T,NULL)); 654 CHKERRQ(MatDestroy(&T)); 655 } 656 CHKERRQ(MatDestroy(&TA)); 657 } 658 659 if (testtt) { /* test with Transpose(Transpose) */ 660 Mat TA, TTA; 661 662 CHKERRQ(MatCreateTranspose(A,&TA)); 663 CHKERRQ(MatCreateTranspose(TA,&TTA)); 664 CHKERRQ(MatMatMult(TTA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X)); 665 CHKERRQ(CheckLocal(B,X,aB,aX)); 666 CHKERRQ(MatMatMultEqual(TTA,B,X,10,&flg)); 667 if (!flg) { 668 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose(Transpose)\n")); 669 CHKERRQ(MatMatMult(TTA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 670 CHKERRQ(MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN)); 671 CHKERRQ(MatView(T,NULL)); 672 CHKERRQ(MatDestroy(&T)); 673 } 674 CHKERRQ(MatDestroy(&TA)); 675 CHKERRQ(MatDestroy(&TTA)); 676 } 677 678 if (testcircular) { /* test circular */ 679 Mat AB; 680 681 CHKERRQ(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AB)); 682 CHKERRQ(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X)); 683 CHKERRQ(CheckLocal(B,X,aB,aX)); 684 if (M == N && N == K) { 685 CHKERRQ(MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B)); 686 } else { 687 CHKERRQ(MatTransposeMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B)); 688 } 689 CHKERRQ(CheckLocal(B,X,aB,aX)); 690 CHKERRQ(MatDestroy(&AB)); 691 } 692 693 /* Test by Pierre Jolivet */ 694 { 695 Mat C,D,D2,AtA; 696 CHKERRQ(MatCreateNormal(A,&AtA)); 697 CHKERRQ(MatDuplicate(X,MAT_DO_NOT_COPY_VALUES,&C)); 698 CHKERRQ(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&D)); 699 CHKERRQ(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&D2)); 700 CHKERRQ(MatSetRandom(B,NULL)); 701 CHKERRQ(MatSetRandom(C,NULL)); 702 CHKERRQ(MatSetRandom(D,NULL)); 703 CHKERRQ(MatSetRandom(D2,NULL)); 704 CHKERRQ(MatProductCreateWithMat(A,B,NULL,C)); 705 CHKERRQ(MatProductSetType(C,MATPRODUCT_AB)); 706 CHKERRQ(MatProductSetFromOptions(C)); 707 CHKERRQ(MatProductSymbolic(C)); 708 CHKERRQ(MatProductCreateWithMat(A,C,NULL,D)); 709 CHKERRQ(MatProductSetType(D, MATPRODUCT_AtB)); 710 CHKERRQ(MatProductSetFromOptions(D)); 711 CHKERRQ(MatProductSymbolic(D)); 712 CHKERRQ(MatProductNumeric(C)); 713 CHKERRQ(MatProductNumeric(D)); 714 CHKERRQ(MatMatMultEqual(AtA,B,D,10,&flg)); 715 if (!flg) { 716 CHKERRQ(MatMatMult(AtA,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T)); 717 CHKERRQ(MatAXPY(T,-1.0,D,SAME_NONZERO_PATTERN)); 718 CHKERRQ(MatView(T,NULL)); 719 CHKERRQ(MatDestroy(&T)); 720 } 721 CHKERRQ(MatDestroy(&C)); 722 CHKERRQ(MatDestroy(&D)); 723 CHKERRQ(MatDestroy(&D2)); 724 CHKERRQ(MatDestroy(&AtA)); 725 } 726 727 CHKERRQ(MatDestroy(&X)); 728 CHKERRQ(MatDestroy(&Bt)); 729 CHKERRQ(MatDestroy(&B)); 730 CHKERRQ(MatDestroy(&A)); 731 CHKERRQ(MatDestroy(&R)); 732 CHKERRQ(MatDestroy(&PtAP)); 733 CHKERRQ(MatDestroy(&RARt)); 734 CHKERRQ(PetscFree(dataX)); 735 CHKERRQ(PetscFree(dataB)); 736 CHKERRQ(PetscFree(dataR)); 737 CHKERRQ(PetscFree(dataBt)); 738 ierr = PetscFinalize(); 739 return ierr; 740 } 741 742 /*TEST 743 744 test: 745 suffix: 1 746 args: -local {{0 1}} -testshellops 747 748 test: 749 output_file: output/ex70_1.out 750 requires: cuda 751 suffix: 1_cuda 752 args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{seqaijcusparse seqaij}} -testshellops {{0 1}} 753 754 test: 755 output_file: output/ex70_1.out 756 nsize: 2 757 suffix: 1_par 758 args: -local {{0 1}} -testmatmatt 0 759 760 test: 761 output_file: output/ex70_1.out 762 requires: cuda 763 nsize: 2 764 suffix: 1_par_cuda 765 args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{mpiaijcusparse mpiaij}} -testnest 0 -testmatmatt 0 -matmatmult_Bbn 3 766 767 test: 768 output_file: output/ex70_1.out 769 suffix: 2 770 nsize: 1 771 args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} 772 773 testset: 774 requires: cuda 775 output_file: output/ex70_1.out 776 nsize: 1 777 args: -M 7 -N 9 -K 2 -local {{0 1}} -testnest 0 -A_mat_type {{seqdensecuda seqdense}} -xgpu {{0 1}} -bgpu {{0 1}} 778 test: 779 requires: !complex 780 suffix: 2_cuda_real 781 test: 782 # complex+single gives a little bigger error in the MatDenseGetColumnVec test 783 requires: complex !single 784 suffix: 2_cuda_complex 785 786 test: 787 output_file: output/ex70_1.out 788 suffix: 2_par 789 nsize: 2 790 args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular -testmatmatt 0 791 792 test: 793 requires: cuda 794 output_file: output/ex70_1.out 795 suffix: 2_par_cuda 796 nsize: 2 797 args: -M 11 -N 9 -K 1 -local {{0 1}} -testcircular 0 -A_mat_type mpiaijcusparse -xgpu -bgpu -testnest 0 -testmatmatt 0 798 799 test: 800 output_file: output/ex70_1.out 801 suffix: 3 802 nsize: {{1 3}} 803 args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type sbaij -symm -testproj 0 -testmatmatt 0 804 805 test: 806 output_file: output/ex70_1.out 807 suffix: 4 808 nsize: 1 809 args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular 810 811 test: 812 output_file: output/ex70_1.out 813 suffix: 5 814 nsize: {{2 4}} 815 args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular -testmatmatt 0 816 817 test: 818 output_file: output/ex70_1.out 819 suffix: 6 820 nsize: 1 821 args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular 822 823 test: 824 output_file: output/ex70_1.out 825 suffix: 7 826 nsize: 1 827 args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type dense -testnest -testcircular 828 829 TEST*/ 830