1 2 /* 3 Defines the basic matrix operations for the SBAIJ (compressed row) 4 matrix storage format. 5 */ 6 #include <../src/mat/impls/baij/seq/baij.h> /*I "petscmat.h" I*/ 7 #include <../src/mat/impls/sbaij/seq/sbaij.h> 8 #include <petscblaslapack.h> 9 10 #include <../src/mat/impls/sbaij/seq/relax.h> 11 #define USESHORT 12 #include <../src/mat/impls/sbaij/seq/relax.h> 13 14 #if defined(PETSC_HAVE_ELEMENTAL) 15 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 16 #endif 17 #if defined(PETSC_HAVE_SCALAPACK) 18 PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 19 #endif 20 PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat,MatType,MatReuse,Mat*); 21 22 /* 23 Checks for missing diagonals 24 */ 25 PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A,PetscBool *missing,PetscInt *dd) 26 { 27 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 28 PetscInt *diag,*ii = a->i,i; 29 30 PetscFunctionBegin; 31 PetscCall(MatMarkDiagonal_SeqSBAIJ(A)); 32 *missing = PETSC_FALSE; 33 if (A->rmap->n > 0 && !ii) { 34 *missing = PETSC_TRUE; 35 if (dd) *dd = 0; 36 PetscCall(PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n")); 37 } else { 38 diag = a->diag; 39 for (i=0; i<a->mbs; i++) { 40 if (diag[i] >= ii[i+1]) { 41 *missing = PETSC_TRUE; 42 if (dd) *dd = i; 43 break; 44 } 45 } 46 } 47 PetscFunctionReturn(0); 48 } 49 50 PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A) 51 { 52 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 53 PetscInt i,j; 54 55 PetscFunctionBegin; 56 if (!a->diag) { 57 PetscCall(PetscMalloc1(a->mbs,&a->diag)); 58 PetscCall(PetscLogObjectMemory((PetscObject)A,a->mbs*sizeof(PetscInt))); 59 a->free_diag = PETSC_TRUE; 60 } 61 for (i=0; i<a->mbs; i++) { 62 a->diag[i] = a->i[i+1]; 63 for (j=a->i[i]; j<a->i[i+1]; j++) { 64 if (a->j[j] == i) { 65 a->diag[i] = j; 66 break; 67 } 68 } 69 } 70 PetscFunctionReturn(0); 71 } 72 73 static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *inia[],const PetscInt *inja[],PetscBool *done) 74 { 75 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 76 PetscInt i,j,n = a->mbs,nz = a->i[n],*tia,*tja,bs = A->rmap->bs,k,l,cnt; 77 PetscInt **ia = (PetscInt**)inia,**ja = (PetscInt**)inja; 78 79 PetscFunctionBegin; 80 *nn = n; 81 if (!ia) PetscFunctionReturn(0); 82 if (symmetric) { 83 PetscCall(MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,PETSC_FALSE,0,0,&tia,&tja)); 84 nz = tia[n]; 85 } else { 86 tia = a->i; tja = a->j; 87 } 88 89 if (!blockcompressed && bs > 1) { 90 (*nn) *= bs; 91 /* malloc & create the natural set of indices */ 92 PetscCall(PetscMalloc1((n+1)*bs,ia)); 93 if (n) { 94 (*ia)[0] = oshift; 95 for (j=1; j<bs; j++) { 96 (*ia)[j] = (tia[1]-tia[0])*bs+(*ia)[j-1]; 97 } 98 } 99 100 for (i=1; i<n; i++) { 101 (*ia)[i*bs] = (tia[i]-tia[i-1])*bs + (*ia)[i*bs-1]; 102 for (j=1; j<bs; j++) { 103 (*ia)[i*bs+j] = (tia[i+1]-tia[i])*bs + (*ia)[i*bs+j-1]; 104 } 105 } 106 if (n) { 107 (*ia)[n*bs] = (tia[n]-tia[n-1])*bs + (*ia)[n*bs-1]; 108 } 109 110 if (inja) { 111 PetscCall(PetscMalloc1(nz*bs*bs,ja)); 112 cnt = 0; 113 for (i=0; i<n; i++) { 114 for (j=0; j<bs; j++) { 115 for (k=tia[i]; k<tia[i+1]; k++) { 116 for (l=0; l<bs; l++) { 117 (*ja)[cnt++] = bs*tja[k] + l; 118 } 119 } 120 } 121 } 122 } 123 124 if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */ 125 PetscCall(PetscFree(tia)); 126 PetscCall(PetscFree(tja)); 127 } 128 } else if (oshift == 1) { 129 if (symmetric) { 130 nz = tia[A->rmap->n/bs]; 131 /* add 1 to i and j indices */ 132 for (i=0; i<A->rmap->n/bs+1; i++) tia[i] = tia[i] + 1; 133 *ia = tia; 134 if (ja) { 135 for (i=0; i<nz; i++) tja[i] = tja[i] + 1; 136 *ja = tja; 137 } 138 } else { 139 nz = a->i[A->rmap->n/bs]; 140 /* malloc space and add 1 to i and j indices */ 141 PetscCall(PetscMalloc1(A->rmap->n/bs+1,ia)); 142 for (i=0; i<A->rmap->n/bs+1; i++) (*ia)[i] = a->i[i] + 1; 143 if (ja) { 144 PetscCall(PetscMalloc1(nz,ja)); 145 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 146 } 147 } 148 } else { 149 *ia = tia; 150 if (ja) *ja = tja; 151 } 152 PetscFunctionReturn(0); 153 } 154 155 static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 156 { 157 PetscFunctionBegin; 158 if (!ia) PetscFunctionReturn(0); 159 if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) { 160 PetscCall(PetscFree(*ia)); 161 if (ja) PetscCall(PetscFree(*ja)); 162 } 163 PetscFunctionReturn(0); 164 } 165 166 PetscErrorCode MatDestroy_SeqSBAIJ(Mat A) 167 { 168 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 169 170 PetscFunctionBegin; 171 #if defined(PETSC_USE_LOG) 172 PetscLogObjectState((PetscObject)A,"Rows=%" PetscInt_FMT ", NZ=%" PetscInt_FMT,A->rmap->N,a->nz); 173 #endif 174 PetscCall(MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i)); 175 if (a->free_diag) PetscCall(PetscFree(a->diag)); 176 PetscCall(ISDestroy(&a->row)); 177 PetscCall(ISDestroy(&a->col)); 178 PetscCall(ISDestroy(&a->icol)); 179 PetscCall(PetscFree(a->idiag)); 180 PetscCall(PetscFree(a->inode.size)); 181 if (a->free_imax_ilen) PetscCall(PetscFree2(a->imax,a->ilen)); 182 PetscCall(PetscFree(a->solve_work)); 183 PetscCall(PetscFree(a->sor_work)); 184 PetscCall(PetscFree(a->solves_work)); 185 PetscCall(PetscFree(a->mult_work)); 186 PetscCall(PetscFree(a->saved_values)); 187 if (a->free_jshort) PetscCall(PetscFree(a->jshort)); 188 PetscCall(PetscFree(a->inew)); 189 PetscCall(MatDestroy(&a->parent)); 190 PetscCall(PetscFree(A->data)); 191 192 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 193 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJGetArray_C",NULL)); 194 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJRestoreArray_C",NULL)); 195 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL)); 196 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL)); 197 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C",NULL)); 198 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C",NULL)); 199 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C",NULL)); 200 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",NULL)); 201 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocationCSR_C",NULL)); 202 #if defined(PETSC_HAVE_ELEMENTAL) 203 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_elemental_C",NULL)); 204 #endif 205 #if defined(PETSC_HAVE_SCALAPACK) 206 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_scalapack_C",NULL)); 207 #endif 208 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL)); 209 PetscFunctionReturn(0); 210 } 211 212 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscBool flg) 213 { 214 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 215 #if defined(PETSC_USE_COMPLEX) 216 PetscInt bs; 217 #endif 218 219 PetscFunctionBegin; 220 #if defined(PETSC_USE_COMPLEX) 221 PetscCall(MatGetBlockSize(A,&bs)); 222 #endif 223 switch (op) { 224 case MAT_ROW_ORIENTED: 225 a->roworiented = flg; 226 break; 227 case MAT_KEEP_NONZERO_PATTERN: 228 a->keepnonzeropattern = flg; 229 break; 230 case MAT_NEW_NONZERO_LOCATIONS: 231 a->nonew = (flg ? 0 : 1); 232 break; 233 case MAT_NEW_NONZERO_LOCATION_ERR: 234 a->nonew = (flg ? -1 : 0); 235 break; 236 case MAT_NEW_NONZERO_ALLOCATION_ERR: 237 a->nonew = (flg ? -2 : 0); 238 break; 239 case MAT_UNUSED_NONZERO_LOCATION_ERR: 240 a->nounused = (flg ? -1 : 0); 241 break; 242 case MAT_FORCE_DIAGONAL_ENTRIES: 243 case MAT_IGNORE_OFF_PROC_ENTRIES: 244 case MAT_USE_HASH_TABLE: 245 case MAT_SORTED_FULL: 246 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 247 break; 248 case MAT_HERMITIAN: 249 #if defined(PETSC_USE_COMPLEX) 250 if (flg) { /* disable transpose ops */ 251 PetscCheck(bs <= 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for Hermitian with block size greater than 1"); 252 A->ops->multtranspose = NULL; 253 A->ops->multtransposeadd = NULL; 254 A->symmetric = PETSC_BOOL3_FALSE; 255 } 256 #endif 257 break; 258 case MAT_SYMMETRIC: 259 case MAT_SPD: 260 #if defined(PETSC_USE_COMPLEX) 261 if (flg) { /* An hermitian and symmetric matrix has zero imaginary part (restore back transpose ops) */ 262 A->ops->multtranspose = A->ops->mult; 263 A->ops->multtransposeadd = A->ops->multadd; 264 } 265 #endif 266 break; 267 /* These options are handled directly by MatSetOption() */ 268 case MAT_STRUCTURALLY_SYMMETRIC: 269 case MAT_SYMMETRY_ETERNAL: 270 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 271 case MAT_STRUCTURE_ONLY: 272 case MAT_SPD_ETERNAL: 273 /* These options are handled directly by MatSetOption() */ 274 break; 275 case MAT_IGNORE_LOWER_TRIANGULAR: 276 a->ignore_ltriangular = flg; 277 break; 278 case MAT_ERROR_LOWER_TRIANGULAR: 279 a->ignore_ltriangular = flg; 280 break; 281 case MAT_GETROW_UPPERTRIANGULAR: 282 a->getrow_utriangular = flg; 283 break; 284 case MAT_SUBMAT_SINGLEIS: 285 break; 286 default: 287 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 288 } 289 PetscFunctionReturn(0); 290 } 291 292 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 293 { 294 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 295 296 PetscFunctionBegin; 297 PetscCheck(!A || a->getrow_utriangular,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE) or MatGetRowUpperTriangular()"); 298 299 /* Get the upper triangular part of the row */ 300 PetscCall(MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a)); 301 PetscFunctionReturn(0); 302 } 303 304 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 305 { 306 PetscFunctionBegin; 307 if (nz) *nz = 0; 308 if (idx) PetscCall(PetscFree(*idx)); 309 if (v) PetscCall(PetscFree(*v)); 310 PetscFunctionReturn(0); 311 } 312 313 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 314 { 315 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 316 317 PetscFunctionBegin; 318 a->getrow_utriangular = PETSC_TRUE; 319 PetscFunctionReturn(0); 320 } 321 322 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 323 { 324 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 325 326 PetscFunctionBegin; 327 a->getrow_utriangular = PETSC_FALSE; 328 PetscFunctionReturn(0); 329 } 330 331 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B) 332 { 333 PetscFunctionBegin; 334 if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A,*B)); 335 if (reuse == MAT_INITIAL_MATRIX) { 336 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,B)); 337 } else if (reuse == MAT_REUSE_MATRIX) { 338 PetscCall(MatCopy(A,*B,SAME_NONZERO_PATTERN)); 339 } 340 PetscFunctionReturn(0); 341 } 342 343 PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 344 { 345 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 346 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 347 PetscViewerFormat format; 348 PetscInt *diag; 349 350 PetscFunctionBegin; 351 PetscCall(PetscViewerGetFormat(viewer,&format)); 352 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 353 PetscCall(PetscViewerASCIIPrintf(viewer," block size is %" PetscInt_FMT "\n",bs)); 354 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 355 Mat aij; 356 const char *matname; 357 358 if (A->factortype && bs>1) { 359 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n")); 360 PetscFunctionReturn(0); 361 } 362 PetscCall(MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij)); 363 PetscCall(PetscObjectGetName((PetscObject)A,&matname)); 364 PetscCall(PetscObjectSetName((PetscObject)aij,matname)); 365 PetscCall(MatView(aij,viewer)); 366 PetscCall(MatDestroy(&aij)); 367 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 368 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 369 for (i=0; i<a->mbs; i++) { 370 for (j=0; j<bs; j++) { 371 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i*bs+j)); 372 for (k=a->i[i]; k<a->i[i+1]; k++) { 373 for (l=0; l<bs; l++) { 374 #if defined(PETSC_USE_COMPLEX) 375 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 376 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k]+l, 377 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 378 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 379 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k]+l, 380 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 381 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 382 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]))); 383 } 384 #else 385 if (a->a[bs2*k + l*bs + j] != 0.0) { 386 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j])); 387 } 388 #endif 389 } 390 } 391 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 392 } 393 } 394 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 395 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 396 PetscFunctionReturn(0); 397 } else { 398 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 399 if (A->factortype) { /* for factored matrix */ 400 PetscCheck(bs<=1,PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet"); 401 402 diag=a->diag; 403 for (i=0; i<a->mbs; i++) { /* for row block i */ 404 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 405 /* diagonal entry */ 406 #if defined(PETSC_USE_COMPLEX) 407 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 408 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),(double)PetscImaginaryPart(1.0/a->a[diag[i]]))); 409 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 410 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),-(double)PetscImaginaryPart(1.0/a->a[diag[i]]))); 411 } else { 412 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]))); 413 } 414 #else 415 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]))); 416 #endif 417 /* off-diagonal entries */ 418 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 419 #if defined(PETSC_USE_COMPLEX) 420 if (PetscImaginaryPart(a->a[k]) > 0.0) { 421 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]))); 422 } else if (PetscImaginaryPart(a->a[k]) < 0.0) { 423 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]))); 424 } else { 425 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]))); 426 } 427 #else 428 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[k],(double)a->a[k])); 429 #endif 430 } 431 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 432 } 433 434 } else { /* for non-factored matrix */ 435 for (i=0; i<a->mbs; i++) { /* for row block i */ 436 for (j=0; j<bs; j++) { /* for row bs*i + j */ 437 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i*bs+j)); 438 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */ 439 for (l=0; l<bs; l++) { /* for column */ 440 #if defined(PETSC_USE_COMPLEX) 441 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 442 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k]+l, 443 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 444 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 445 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k]+l, 446 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 447 } else { 448 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]))); 449 } 450 #else 451 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j])); 452 #endif 453 } 454 } 455 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 456 } 457 } 458 } 459 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 460 } 461 PetscCall(PetscViewerFlush(viewer)); 462 PetscFunctionReturn(0); 463 } 464 465 #include <petscdraw.h> 466 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 467 { 468 Mat A = (Mat) Aa; 469 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 470 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 471 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 472 MatScalar *aa; 473 PetscViewer viewer; 474 475 PetscFunctionBegin; 476 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 477 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 478 479 /* loop over matrix elements drawing boxes */ 480 481 PetscDrawCollectiveBegin(draw); 482 PetscCall(PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric")); 483 /* Blue for negative, Cyan for zero and Red for positive */ 484 color = PETSC_DRAW_BLUE; 485 for (i=0,row=0; i<mbs; i++,row+=bs) { 486 for (j=a->i[i]; j<a->i[i+1]; j++) { 487 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 488 x_l = a->j[j]*bs; x_r = x_l + 1.0; 489 aa = a->a + j*bs2; 490 for (k=0; k<bs; k++) { 491 for (l=0; l<bs; l++) { 492 if (PetscRealPart(*aa++) >= 0.) continue; 493 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 494 } 495 } 496 } 497 } 498 color = PETSC_DRAW_CYAN; 499 for (i=0,row=0; i<mbs; i++,row+=bs) { 500 for (j=a->i[i]; j<a->i[i+1]; j++) { 501 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 502 x_l = a->j[j]*bs; x_r = x_l + 1.0; 503 aa = a->a + j*bs2; 504 for (k=0; k<bs; k++) { 505 for (l=0; l<bs; l++) { 506 if (PetscRealPart(*aa++) != 0.) continue; 507 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 508 } 509 } 510 } 511 } 512 color = PETSC_DRAW_RED; 513 for (i=0,row=0; i<mbs; i++,row+=bs) { 514 for (j=a->i[i]; j<a->i[i+1]; j++) { 515 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 516 x_l = a->j[j]*bs; x_r = x_l + 1.0; 517 aa = a->a + j*bs2; 518 for (k=0; k<bs; k++) { 519 for (l=0; l<bs; l++) { 520 if (PetscRealPart(*aa++) <= 0.) continue; 521 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 522 } 523 } 524 } 525 } 526 PetscDrawCollectiveEnd(draw); 527 PetscFunctionReturn(0); 528 } 529 530 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 531 { 532 PetscReal xl,yl,xr,yr,w,h; 533 PetscDraw draw; 534 PetscBool isnull; 535 536 PetscFunctionBegin; 537 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 538 PetscCall(PetscDrawIsNull(draw,&isnull)); 539 if (isnull) PetscFunctionReturn(0); 540 541 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 542 xr += w; yr += h; xl = -w; yl = -h; 543 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 544 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 545 PetscCall(PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A)); 546 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 547 PetscCall(PetscDrawSave(draw)); 548 PetscFunctionReturn(0); 549 } 550 551 /* Used for both MPIBAIJ and MPISBAIJ matrices */ 552 #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary 553 554 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 555 { 556 PetscBool iascii,isbinary,isdraw; 557 558 PetscFunctionBegin; 559 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 560 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 561 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 562 if (iascii) { 563 PetscCall(MatView_SeqSBAIJ_ASCII(A,viewer)); 564 } else if (isbinary) { 565 PetscCall(MatView_SeqSBAIJ_Binary(A,viewer)); 566 } else if (isdraw) { 567 PetscCall(MatView_SeqSBAIJ_Draw(A,viewer)); 568 } else { 569 Mat B; 570 const char *matname; 571 PetscCall(MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B)); 572 PetscCall(PetscObjectGetName((PetscObject)A,&matname)); 573 PetscCall(PetscObjectSetName((PetscObject)B,matname)); 574 PetscCall(MatView(B,viewer)); 575 PetscCall(MatDestroy(&B)); 576 } 577 PetscFunctionReturn(0); 578 } 579 580 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 581 { 582 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 583 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 584 PetscInt *ai = a->i,*ailen = a->ilen; 585 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 586 MatScalar *ap,*aa = a->a; 587 588 PetscFunctionBegin; 589 for (k=0; k<m; k++) { /* loop over rows */ 590 row = im[k]; brow = row/bs; 591 if (row < 0) {v += n; continue;} /* negative row */ 592 PetscCheck(row < A->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->N-1); 593 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 594 nrow = ailen[brow]; 595 for (l=0; l<n; l++) { /* loop over columns */ 596 if (in[l] < 0) {v++; continue;} /* negative column */ 597 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 598 col = in[l]; 599 bcol = col/bs; 600 cidx = col%bs; 601 ridx = row%bs; 602 high = nrow; 603 low = 0; /* assume unsorted */ 604 while (high-low > 5) { 605 t = (low+high)/2; 606 if (rp[t] > bcol) high = t; 607 else low = t; 608 } 609 for (i=low; i<high; i++) { 610 if (rp[i] > bcol) break; 611 if (rp[i] == bcol) { 612 *v++ = ap[bs2*i+bs*cidx+ridx]; 613 goto finished; 614 } 615 } 616 *v++ = 0.0; 617 finished:; 618 } 619 } 620 PetscFunctionReturn(0); 621 } 622 623 PetscErrorCode MatPermute_SeqSBAIJ(Mat A,IS rowp,IS colp,Mat *B) 624 { 625 Mat C; 626 627 PetscFunctionBegin; 628 PetscCall(MatConvert(A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&C)); 629 PetscCall(MatPermute(C,rowp,colp,B)); 630 PetscCall(MatDestroy(&C)); 631 if (rowp == colp) { 632 PetscCall(MatConvert(*B,MATSEQSBAIJ,MAT_INPLACE_MATRIX,B)); 633 } 634 PetscFunctionReturn(0); 635 } 636 637 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 638 { 639 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 640 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 641 PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen; 642 PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 643 PetscBool roworiented=a->roworiented; 644 const PetscScalar *value = v; 645 MatScalar *ap,*aa = a->a,*bap; 646 647 PetscFunctionBegin; 648 if (roworiented) stepval = (n-1)*bs; 649 else stepval = (m-1)*bs; 650 651 for (k=0; k<m; k++) { /* loop over added rows */ 652 row = im[k]; 653 if (row < 0) continue; 654 PetscCheck(row < a->mbs,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index row too large %" PetscInt_FMT " max %" PetscInt_FMT,row,a->mbs-1); 655 rp = aj + ai[row]; 656 ap = aa + bs2*ai[row]; 657 rmax = imax[row]; 658 nrow = ailen[row]; 659 low = 0; 660 high = nrow; 661 for (l=0; l<n; l++) { /* loop over added columns */ 662 if (in[l] < 0) continue; 663 col = in[l]; 664 PetscCheck(col < a->nbs,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index column too large %" PetscInt_FMT " max %" PetscInt_FMT,col,a->nbs-1); 665 if (col < row) { 666 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */ 667 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 668 } 669 if (roworiented) value = v + k*(stepval+bs)*bs + l*bs; 670 else value = v + l*(stepval+bs)*bs + k*bs; 671 672 if (col <= lastcol) low = 0; 673 else high = nrow; 674 675 lastcol = col; 676 while (high-low > 7) { 677 t = (low+high)/2; 678 if (rp[t] > col) high = t; 679 else low = t; 680 } 681 for (i=low; i<high; i++) { 682 if (rp[i] > col) break; 683 if (rp[i] == col) { 684 bap = ap + bs2*i; 685 if (roworiented) { 686 if (is == ADD_VALUES) { 687 for (ii=0; ii<bs; ii++,value+=stepval) { 688 for (jj=ii; jj<bs2; jj+=bs) { 689 bap[jj] += *value++; 690 } 691 } 692 } else { 693 for (ii=0; ii<bs; ii++,value+=stepval) { 694 for (jj=ii; jj<bs2; jj+=bs) { 695 bap[jj] = *value++; 696 } 697 } 698 } 699 } else { 700 if (is == ADD_VALUES) { 701 for (ii=0; ii<bs; ii++,value+=stepval) { 702 for (jj=0; jj<bs; jj++) { 703 *bap++ += *value++; 704 } 705 } 706 } else { 707 for (ii=0; ii<bs; ii++,value+=stepval) { 708 for (jj=0; jj<bs; jj++) { 709 *bap++ = *value++; 710 } 711 } 712 } 713 } 714 goto noinsert2; 715 } 716 } 717 if (nonew == 1) goto noinsert2; 718 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new block index nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col); 719 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 720 N = nrow++ - 1; high++; 721 /* shift up all the later entries in this row */ 722 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 723 PetscCall(PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1))); 724 PetscCall(PetscArrayzero(ap+bs2*i,bs2)); 725 rp[i] = col; 726 bap = ap + bs2*i; 727 if (roworiented) { 728 for (ii=0; ii<bs; ii++,value+=stepval) { 729 for (jj=ii; jj<bs2; jj+=bs) { 730 bap[jj] = *value++; 731 } 732 } 733 } else { 734 for (ii=0; ii<bs; ii++,value+=stepval) { 735 for (jj=0; jj<bs; jj++) { 736 *bap++ = *value++; 737 } 738 } 739 } 740 noinsert2:; 741 low = i; 742 } 743 ailen[row] = nrow; 744 } 745 PetscFunctionReturn(0); 746 } 747 748 /* 749 This is not yet used 750 */ 751 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 752 { 753 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 754 const PetscInt *ai = a->i, *aj = a->j,*cols; 755 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 756 PetscBool flag; 757 758 PetscFunctionBegin; 759 PetscCall(PetscMalloc1(m,&ns)); 760 while (i < m) { 761 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 762 /* Limits the number of elements in a node to 'a->inode.limit' */ 763 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 764 nzy = ai[j+1] - ai[j]; 765 if (nzy != (nzx - j + i)) break; 766 PetscCall(PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag)); 767 if (!flag) break; 768 } 769 ns[node_count++] = blk_size; 770 771 i = j; 772 } 773 if (!a->inode.size && m && node_count > .9*m) { 774 PetscCall(PetscFree(ns)); 775 PetscCall(PetscInfo(A,"Found %" PetscInt_FMT " nodes out of %" PetscInt_FMT " rows. Not using Inode routines\n",node_count,m)); 776 } else { 777 a->inode.node_count = node_count; 778 779 PetscCall(PetscMalloc1(node_count,&a->inode.size)); 780 PetscCall(PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt))); 781 PetscCall(PetscArraycpy(a->inode.size,ns,node_count)); 782 PetscCall(PetscFree(ns)); 783 PetscCall(PetscInfo(A,"Found %" PetscInt_FMT " nodes of %" PetscInt_FMT ". Limit used: %" PetscInt_FMT ". Using Inode routines\n",node_count,m,a->inode.limit)); 784 785 /* count collections of adjacent columns in each inode */ 786 row = 0; 787 cnt = 0; 788 for (i=0; i<node_count; i++) { 789 cols = aj + ai[row] + a->inode.size[i]; 790 nz = ai[row+1] - ai[row] - a->inode.size[i]; 791 for (j=1; j<nz; j++) { 792 if (cols[j] != cols[j-1]+1) cnt++; 793 } 794 cnt++; 795 row += a->inode.size[i]; 796 } 797 PetscCall(PetscMalloc1(2*cnt,&counts)); 798 cnt = 0; 799 row = 0; 800 for (i=0; i<node_count; i++) { 801 cols = aj + ai[row] + a->inode.size[i]; 802 counts[2*cnt] = cols[0]; 803 nz = ai[row+1] - ai[row] - a->inode.size[i]; 804 cnt2 = 1; 805 for (j=1; j<nz; j++) { 806 if (cols[j] != cols[j-1]+1) { 807 counts[2*(cnt++)+1] = cnt2; 808 counts[2*cnt] = cols[j]; 809 cnt2 = 1; 810 } else cnt2++; 811 } 812 counts[2*(cnt++)+1] = cnt2; 813 row += a->inode.size[i]; 814 } 815 PetscCall(PetscIntView(2*cnt,counts,NULL)); 816 } 817 PetscFunctionReturn(0); 818 } 819 820 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 821 { 822 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 823 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 824 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 825 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 826 MatScalar *aa = a->a,*ap; 827 828 PetscFunctionBegin; 829 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 830 831 if (m) rmax = ailen[0]; 832 for (i=1; i<mbs; i++) { 833 /* move each row back by the amount of empty slots (fshift) before it*/ 834 fshift += imax[i-1] - ailen[i-1]; 835 rmax = PetscMax(rmax,ailen[i]); 836 if (fshift) { 837 ip = aj + ai[i]; 838 ap = aa + bs2*ai[i]; 839 N = ailen[i]; 840 PetscCall(PetscArraymove(ip-fshift,ip,N)); 841 PetscCall(PetscArraymove(ap-bs2*fshift,ap,bs2*N)); 842 } 843 ai[i] = ai[i-1] + ailen[i-1]; 844 } 845 if (mbs) { 846 fshift += imax[mbs-1] - ailen[mbs-1]; 847 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 848 } 849 /* reset ilen and imax for each row */ 850 for (i=0; i<mbs; i++) { 851 ailen[i] = imax[i] = ai[i+1] - ai[i]; 852 } 853 a->nz = ai[mbs]; 854 855 /* diagonals may have moved, reset it */ 856 if (a->diag) PetscCall(PetscArraycpy(a->diag,ai,mbs)); 857 PetscCheck(!fshift || a->nounused != -1,PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT " block size %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, A->rmap->bs, fshift*bs2); 858 859 PetscCall(PetscInfo(A,"Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT ", block size %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded, %" PetscInt_FMT " used\n",m,A->rmap->N,A->rmap->bs,fshift*bs2,a->nz*bs2)); 860 PetscCall(PetscInfo(A,"Number of mallocs during MatSetValues is %" PetscInt_FMT "\n",a->reallocs)); 861 PetscCall(PetscInfo(A,"Most nonzeros blocks in any row is %" PetscInt_FMT "\n",rmax)); 862 863 A->info.mallocs += a->reallocs; 864 a->reallocs = 0; 865 A->info.nz_unneeded = (PetscReal)fshift*bs2; 866 a->idiagvalid = PETSC_FALSE; 867 a->rmax = rmax; 868 869 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 870 if (a->jshort && a->free_jshort) { 871 /* when matrix data structure is changed, previous jshort must be replaced */ 872 PetscCall(PetscFree(a->jshort)); 873 } 874 PetscCall(PetscMalloc1(a->i[A->rmap->n],&a->jshort)); 875 PetscCall(PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short))); 876 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 877 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 878 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 879 a->free_jshort = PETSC_TRUE; 880 } 881 PetscFunctionReturn(0); 882 } 883 884 /* 885 This function returns an array of flags which indicate the locations of contiguous 886 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 887 then the resulting sizes = [3,1,1,3,1] corresponding to sets [(0,1,2),(3),(5),(6,7,8),(9)] 888 Assume: sizes should be long enough to hold all the values. 889 */ 890 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 891 { 892 PetscInt i,j,k,row; 893 PetscBool flg; 894 895 PetscFunctionBegin; 896 for (i=0,j=0; i<n; j++) { 897 row = idx[i]; 898 if (row%bs!=0) { /* Not the beginning of a block */ 899 sizes[j] = 1; 900 i++; 901 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 902 sizes[j] = 1; /* Also makes sure at least 'bs' values exist for next else */ 903 i++; 904 } else { /* Beginning of the block, so check if the complete block exists */ 905 flg = PETSC_TRUE; 906 for (k=1; k<bs; k++) { 907 if (row+k != idx[i+k]) { /* break in the block */ 908 flg = PETSC_FALSE; 909 break; 910 } 911 } 912 if (flg) { /* No break in the bs */ 913 sizes[j] = bs; 914 i += bs; 915 } else { 916 sizes[j] = 1; 917 i++; 918 } 919 } 920 } 921 *bs_max = j; 922 PetscFunctionReturn(0); 923 } 924 925 /* Only add/insert a(i,j) with i<=j (blocks). 926 Any a(i,j) with i>j input by user is ingored. 927 */ 928 929 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 930 { 931 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 932 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 933 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 934 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 935 PetscInt ridx,cidx,bs2=a->bs2; 936 MatScalar *ap,value,*aa=a->a,*bap; 937 938 PetscFunctionBegin; 939 for (k=0; k<m; k++) { /* loop over added rows */ 940 row = im[k]; /* row number */ 941 brow = row/bs; /* block row number */ 942 if (row < 0) continue; 943 PetscCheck(row < A->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->N-1); 944 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 945 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 946 rmax = imax[brow]; /* maximum space allocated for this row */ 947 nrow = ailen[brow]; /* actual length of this row */ 948 low = 0; 949 high = nrow; 950 for (l=0; l<n; l++) { /* loop over added columns */ 951 if (in[l] < 0) continue; 952 PetscCheck(in[l] < A->cmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->N-1); 953 col = in[l]; 954 bcol = col/bs; /* block col number */ 955 956 if (brow > bcol) { 957 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */ 958 else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 959 } 960 961 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 962 if ((brow==bcol && ridx<=cidx) || (brow<bcol)) { 963 /* element value a(k,l) */ 964 if (roworiented) value = v[l + k*n]; 965 else value = v[k + l*m]; 966 967 /* move pointer bap to a(k,l) quickly and add/insert value */ 968 if (col <= lastcol) low = 0; 969 else high = nrow; 970 971 lastcol = col; 972 while (high-low > 7) { 973 t = (low+high)/2; 974 if (rp[t] > bcol) high = t; 975 else low = t; 976 } 977 for (i=low; i<high; i++) { 978 if (rp[i] > bcol) break; 979 if (rp[i] == bcol) { 980 bap = ap + bs2*i + bs*cidx + ridx; 981 if (is == ADD_VALUES) *bap += value; 982 else *bap = value; 983 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 984 if (brow == bcol && ridx < cidx) { 985 bap = ap + bs2*i + bs*ridx + cidx; 986 if (is == ADD_VALUES) *bap += value; 987 else *bap = value; 988 } 989 goto noinsert1; 990 } 991 } 992 993 if (nonew == 1) goto noinsert1; 994 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col); 995 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 996 997 N = nrow++ - 1; high++; 998 /* shift up all the later entries in this row */ 999 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 1000 PetscCall(PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1))); 1001 PetscCall(PetscArrayzero(ap+bs2*i,bs2)); 1002 rp[i] = bcol; 1003 ap[bs2*i + bs*cidx + ridx] = value; 1004 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1005 if (brow == bcol && ridx < cidx) { 1006 ap[bs2*i + bs*ridx + cidx] = value; 1007 } 1008 A->nonzerostate++; 1009 noinsert1:; 1010 low = i; 1011 } 1012 } /* end of loop over added columns */ 1013 ailen[brow] = nrow; 1014 } /* end of loop over added rows */ 1015 PetscFunctionReturn(0); 1016 } 1017 1018 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1019 { 1020 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1021 Mat outA; 1022 PetscBool row_identity; 1023 1024 PetscFunctionBegin; 1025 PetscCheck(info->levels == 0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1026 PetscCall(ISIdentity(row,&row_identity)); 1027 PetscCheck(row_identity,PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1028 PetscCheck(inA->rmap->bs == 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix block size %" PetscInt_FMT " is not supported",inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */ 1029 1030 outA = inA; 1031 inA->factortype = MAT_FACTOR_ICC; 1032 PetscCall(PetscFree(inA->solvertype)); 1033 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype)); 1034 1035 PetscCall(MatMarkDiagonal_SeqSBAIJ(inA)); 1036 PetscCall(MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity)); 1037 1038 PetscCall(PetscObjectReference((PetscObject)row)); 1039 PetscCall(ISDestroy(&a->row)); 1040 a->row = row; 1041 PetscCall(PetscObjectReference((PetscObject)row)); 1042 PetscCall(ISDestroy(&a->col)); 1043 a->col = row; 1044 1045 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1046 if (a->icol) PetscCall(ISInvertPermutation(row,PETSC_DECIDE, &a->icol)); 1047 PetscCall(PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol)); 1048 1049 if (!a->solve_work) { 1050 PetscCall(PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work)); 1051 PetscCall(PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar))); 1052 } 1053 1054 PetscCall(MatCholeskyFactorNumeric(outA,inA,info)); 1055 PetscFunctionReturn(0); 1056 } 1057 1058 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1059 { 1060 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data; 1061 PetscInt i,nz,n; 1062 1063 PetscFunctionBegin; 1064 nz = baij->maxnz; 1065 n = mat->cmap->n; 1066 for (i=0; i<nz; i++) baij->j[i] = indices[i]; 1067 1068 baij->nz = nz; 1069 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i]; 1070 1071 PetscCall(MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 1072 PetscFunctionReturn(0); 1073 } 1074 1075 /*@ 1076 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1077 in the matrix. 1078 1079 Input Parameters: 1080 + mat - the SeqSBAIJ matrix 1081 - indices - the column indices 1082 1083 Level: advanced 1084 1085 Notes: 1086 This can be called if you have precomputed the nonzero structure of the 1087 matrix and want to provide it to the matrix object to improve the performance 1088 of the MatSetValues() operation. 1089 1090 You MUST have set the correct numbers of nonzeros per row in the call to 1091 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1092 1093 MUST be called before any calls to MatSetValues() 1094 1095 .seealso: `MatCreateSeqSBAIJ` 1096 @*/ 1097 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1098 { 1099 PetscFunctionBegin; 1100 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1101 PetscValidIntPointer(indices,2); 1102 PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices)); 1103 PetscFunctionReturn(0); 1104 } 1105 1106 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1107 { 1108 PetscBool isbaij; 1109 1110 PetscFunctionBegin; 1111 PetscCall(PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"")); 1112 PetscCheck(isbaij,PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name); 1113 /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */ 1114 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 1115 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1116 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1117 1118 PetscCheck(a->i[a->mbs] == b->i[b->mbs],PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1119 PetscCheck(a->mbs == b->mbs,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different"); 1120 PetscCheck(a->bs2 == b->bs2,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size"); 1121 PetscCall(PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs])); 1122 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 1123 } else { 1124 PetscCall(MatGetRowUpperTriangular(A)); 1125 PetscCall(MatCopy_Basic(A,B,str)); 1126 PetscCall(MatRestoreRowUpperTriangular(A)); 1127 } 1128 PetscFunctionReturn(0); 1129 } 1130 1131 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) 1132 { 1133 PetscFunctionBegin; 1134 PetscCall(MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL)); 1135 PetscFunctionReturn(0); 1136 } 1137 1138 static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1139 { 1140 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1141 1142 PetscFunctionBegin; 1143 *array = a->a; 1144 PetscFunctionReturn(0); 1145 } 1146 1147 static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1148 { 1149 PetscFunctionBegin; 1150 *array = NULL; 1151 PetscFunctionReturn(0); 1152 } 1153 1154 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz) 1155 { 1156 PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs; 1157 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data; 1158 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data; 1159 1160 PetscFunctionBegin; 1161 /* Set the number of nonzeros in the new matrix */ 1162 PetscCall(MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz)); 1163 PetscFunctionReturn(0); 1164 } 1165 1166 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1167 { 1168 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data; 1169 PetscInt bs=Y->rmap->bs,bs2=bs*bs; 1170 PetscBLASInt one = 1; 1171 1172 PetscFunctionBegin; 1173 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 1174 PetscBool e = x->nz == y->nz && x->mbs == y->mbs ? PETSC_TRUE : PETSC_FALSE; 1175 if (e) { 1176 PetscCall(PetscArraycmp(x->i,y->i,x->mbs+1,&e)); 1177 if (e) { 1178 PetscCall(PetscArraycmp(x->j,y->j,x->i[x->mbs],&e)); 1179 if (e) str = SAME_NONZERO_PATTERN; 1180 } 1181 } 1182 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatStructure is not SAME_NONZERO_PATTERN"); 1183 } 1184 if (str == SAME_NONZERO_PATTERN) { 1185 PetscScalar alpha = a; 1186 PetscBLASInt bnz; 1187 PetscCall(PetscBLASIntCast(x->nz*bs2,&bnz)); 1188 PetscCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 1189 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 1190 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1191 PetscCall(MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE)); 1192 PetscCall(MatAXPY_Basic(Y,a,X,str)); 1193 PetscCall(MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE)); 1194 } else { 1195 Mat B; 1196 PetscInt *nnz; 1197 PetscCheck(bs == X->rmap->bs,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size"); 1198 PetscCall(MatGetRowUpperTriangular(X)); 1199 PetscCall(MatGetRowUpperTriangular(Y)); 1200 PetscCall(PetscMalloc1(Y->rmap->N,&nnz)); 1201 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),&B)); 1202 PetscCall(PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name)); 1203 PetscCall(MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N)); 1204 PetscCall(MatSetBlockSizesFromMats(B,Y,Y)); 1205 PetscCall(MatSetType(B,((PetscObject)Y)->type_name)); 1206 PetscCall(MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz)); 1207 PetscCall(MatSeqSBAIJSetPreallocation(B,bs,0,nnz)); 1208 1209 PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 1210 1211 PetscCall(MatHeaderMerge(Y,&B)); 1212 PetscCall(PetscFree(nnz)); 1213 PetscCall(MatRestoreRowUpperTriangular(X)); 1214 PetscCall(MatRestoreRowUpperTriangular(Y)); 1215 } 1216 PetscFunctionReturn(0); 1217 } 1218 1219 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1220 { 1221 PetscFunctionBegin; 1222 *flg = PETSC_TRUE; 1223 PetscFunctionReturn(0); 1224 } 1225 1226 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1227 { 1228 PetscFunctionBegin; 1229 *flg = PETSC_TRUE; 1230 PetscFunctionReturn(0); 1231 } 1232 1233 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1234 { 1235 PetscFunctionBegin; 1236 *flg = PETSC_FALSE; 1237 PetscFunctionReturn(0); 1238 } 1239 1240 PetscErrorCode MatConjugate_SeqSBAIJ(Mat A) 1241 { 1242 #if defined(PETSC_USE_COMPLEX) 1243 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1244 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1245 MatScalar *aa = a->a; 1246 1247 PetscFunctionBegin; 1248 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1249 #else 1250 PetscFunctionBegin; 1251 #endif 1252 PetscFunctionReturn(0); 1253 } 1254 1255 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1256 { 1257 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1258 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1259 MatScalar *aa = a->a; 1260 1261 PetscFunctionBegin; 1262 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1263 PetscFunctionReturn(0); 1264 } 1265 1266 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1267 { 1268 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1269 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1270 MatScalar *aa = a->a; 1271 1272 PetscFunctionBegin; 1273 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1274 PetscFunctionReturn(0); 1275 } 1276 1277 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1278 { 1279 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1280 PetscInt i,j,k,count; 1281 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col; 1282 PetscScalar zero = 0.0; 1283 MatScalar *aa; 1284 const PetscScalar *xx; 1285 PetscScalar *bb; 1286 PetscBool *zeroed,vecs = PETSC_FALSE; 1287 1288 PetscFunctionBegin; 1289 /* fix right hand side if needed */ 1290 if (x && b) { 1291 PetscCall(VecGetArrayRead(x,&xx)); 1292 PetscCall(VecGetArray(b,&bb)); 1293 vecs = PETSC_TRUE; 1294 } 1295 1296 /* zero the columns */ 1297 PetscCall(PetscCalloc1(A->rmap->n,&zeroed)); 1298 for (i=0; i<is_n; i++) { 1299 PetscCheck(is_idx[i] >= 0 && is_idx[i] < A->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range",is_idx[i]); 1300 zeroed[is_idx[i]] = PETSC_TRUE; 1301 } 1302 if (vecs) { 1303 for (i=0; i<A->rmap->N; i++) { 1304 row = i/bs; 1305 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1306 for (k=0; k<bs; k++) { 1307 col = bs*baij->j[j] + k; 1308 if (col <= i) continue; 1309 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1310 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1311 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1312 } 1313 } 1314 } 1315 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1316 } 1317 1318 for (i=0; i<A->rmap->N; i++) { 1319 if (!zeroed[i]) { 1320 row = i/bs; 1321 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1322 for (k=0; k<bs; k++) { 1323 col = bs*baij->j[j] + k; 1324 if (zeroed[col]) { 1325 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1326 aa[0] = 0.0; 1327 } 1328 } 1329 } 1330 } 1331 } 1332 PetscCall(PetscFree(zeroed)); 1333 if (vecs) { 1334 PetscCall(VecRestoreArrayRead(x,&xx)); 1335 PetscCall(VecRestoreArray(b,&bb)); 1336 } 1337 1338 /* zero the rows */ 1339 for (i=0; i<is_n; i++) { 1340 row = is_idx[i]; 1341 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1342 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1343 for (k=0; k<count; k++) { 1344 aa[0] = zero; 1345 aa += bs; 1346 } 1347 if (diag != 0.0) PetscUseTypeMethod(A,setvalues ,1,&row,1,&row,&diag,INSERT_VALUES); 1348 } 1349 PetscCall(MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY)); 1350 PetscFunctionReturn(0); 1351 } 1352 1353 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a) 1354 { 1355 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data; 1356 1357 PetscFunctionBegin; 1358 if (!Y->preallocated || !aij->nz) { 1359 PetscCall(MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL)); 1360 } 1361 PetscCall(MatShift_Basic(Y,a)); 1362 PetscFunctionReturn(0); 1363 } 1364 1365 /* -------------------------------------------------------------------*/ 1366 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1367 MatGetRow_SeqSBAIJ, 1368 MatRestoreRow_SeqSBAIJ, 1369 MatMult_SeqSBAIJ_N, 1370 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1371 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1372 MatMultAdd_SeqSBAIJ_N, 1373 NULL, 1374 NULL, 1375 NULL, 1376 /* 10*/ NULL, 1377 NULL, 1378 MatCholeskyFactor_SeqSBAIJ, 1379 MatSOR_SeqSBAIJ, 1380 MatTranspose_SeqSBAIJ, 1381 /* 15*/ MatGetInfo_SeqSBAIJ, 1382 MatEqual_SeqSBAIJ, 1383 MatGetDiagonal_SeqSBAIJ, 1384 MatDiagonalScale_SeqSBAIJ, 1385 MatNorm_SeqSBAIJ, 1386 /* 20*/ NULL, 1387 MatAssemblyEnd_SeqSBAIJ, 1388 MatSetOption_SeqSBAIJ, 1389 MatZeroEntries_SeqSBAIJ, 1390 /* 24*/ NULL, 1391 NULL, 1392 NULL, 1393 NULL, 1394 NULL, 1395 /* 29*/ MatSetUp_SeqSBAIJ, 1396 NULL, 1397 NULL, 1398 NULL, 1399 NULL, 1400 /* 34*/ MatDuplicate_SeqSBAIJ, 1401 NULL, 1402 NULL, 1403 NULL, 1404 MatICCFactor_SeqSBAIJ, 1405 /* 39*/ MatAXPY_SeqSBAIJ, 1406 MatCreateSubMatrices_SeqSBAIJ, 1407 MatIncreaseOverlap_SeqSBAIJ, 1408 MatGetValues_SeqSBAIJ, 1409 MatCopy_SeqSBAIJ, 1410 /* 44*/ NULL, 1411 MatScale_SeqSBAIJ, 1412 MatShift_SeqSBAIJ, 1413 NULL, 1414 MatZeroRowsColumns_SeqSBAIJ, 1415 /* 49*/ NULL, 1416 MatGetRowIJ_SeqSBAIJ, 1417 MatRestoreRowIJ_SeqSBAIJ, 1418 NULL, 1419 NULL, 1420 /* 54*/ NULL, 1421 NULL, 1422 NULL, 1423 MatPermute_SeqSBAIJ, 1424 MatSetValuesBlocked_SeqSBAIJ, 1425 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1426 NULL, 1427 NULL, 1428 NULL, 1429 NULL, 1430 /* 64*/ NULL, 1431 NULL, 1432 NULL, 1433 NULL, 1434 NULL, 1435 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1436 NULL, 1437 MatConvert_MPISBAIJ_Basic, 1438 NULL, 1439 NULL, 1440 /* 74*/ NULL, 1441 NULL, 1442 NULL, 1443 NULL, 1444 NULL, 1445 /* 79*/ NULL, 1446 NULL, 1447 NULL, 1448 MatGetInertia_SeqSBAIJ, 1449 MatLoad_SeqSBAIJ, 1450 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1451 MatIsHermitian_SeqSBAIJ, 1452 MatIsStructurallySymmetric_SeqSBAIJ, 1453 NULL, 1454 NULL, 1455 /* 89*/ NULL, 1456 NULL, 1457 NULL, 1458 NULL, 1459 NULL, 1460 /* 94*/ NULL, 1461 NULL, 1462 NULL, 1463 NULL, 1464 NULL, 1465 /* 99*/ NULL, 1466 NULL, 1467 NULL, 1468 MatConjugate_SeqSBAIJ, 1469 NULL, 1470 /*104*/ NULL, 1471 MatRealPart_SeqSBAIJ, 1472 MatImaginaryPart_SeqSBAIJ, 1473 MatGetRowUpperTriangular_SeqSBAIJ, 1474 MatRestoreRowUpperTriangular_SeqSBAIJ, 1475 /*109*/ NULL, 1476 NULL, 1477 NULL, 1478 NULL, 1479 MatMissingDiagonal_SeqSBAIJ, 1480 /*114*/ NULL, 1481 NULL, 1482 NULL, 1483 NULL, 1484 NULL, 1485 /*119*/ NULL, 1486 NULL, 1487 NULL, 1488 NULL, 1489 NULL, 1490 /*124*/ NULL, 1491 NULL, 1492 NULL, 1493 NULL, 1494 NULL, 1495 /*129*/ NULL, 1496 NULL, 1497 NULL, 1498 NULL, 1499 NULL, 1500 /*134*/ NULL, 1501 NULL, 1502 NULL, 1503 NULL, 1504 NULL, 1505 /*139*/ MatSetBlockSizes_Default, 1506 NULL, 1507 NULL, 1508 NULL, 1509 NULL, 1510 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ, 1511 NULL, 1512 NULL, 1513 NULL, 1514 NULL, 1515 NULL, 1516 /*150*/ NULL 1517 }; 1518 1519 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1520 { 1521 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1522 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1523 1524 PetscFunctionBegin; 1525 PetscCheck(aij->nonew == 1,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1526 1527 /* allocate space for values if not already there */ 1528 if (!aij->saved_values) { 1529 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 1530 } 1531 1532 /* copy values over */ 1533 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 1534 PetscFunctionReturn(0); 1535 } 1536 1537 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1538 { 1539 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1540 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1541 1542 PetscFunctionBegin; 1543 PetscCheck(aij->nonew == 1,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1544 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1545 1546 /* copy values over */ 1547 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 1548 PetscFunctionReturn(0); 1549 } 1550 1551 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1552 { 1553 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1554 PetscInt i,mbs,nbs,bs2; 1555 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1556 1557 PetscFunctionBegin; 1558 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1559 1560 PetscCall(MatSetBlockSize(B,PetscAbs(bs))); 1561 PetscCall(PetscLayoutSetUp(B->rmap)); 1562 PetscCall(PetscLayoutSetUp(B->cmap)); 1563 PetscCheck(B->rmap->N <= B->cmap->N,PETSC_COMM_SELF,PETSC_ERR_SUP,"SEQSBAIJ matrix cannot have more rows %" PetscInt_FMT " than columns %" PetscInt_FMT,B->rmap->N,B->cmap->N); 1564 PetscCall(PetscLayoutGetBlockSize(B->rmap,&bs)); 1565 1566 B->preallocated = PETSC_TRUE; 1567 1568 mbs = B->rmap->N/bs; 1569 nbs = B->cmap->n/bs; 1570 bs2 = bs*bs; 1571 1572 PetscCheck(mbs*bs == B->rmap->N && nbs*bs == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1573 1574 if (nz == MAT_SKIP_ALLOCATION) { 1575 skipallocation = PETSC_TRUE; 1576 nz = 0; 1577 } 1578 1579 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1580 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 1581 if (nnz) { 1582 for (i=0; i<mbs; i++) { 1583 PetscCheck(nnz[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT,i,nnz[i]); 1584 PetscCheck(nnz[i] <= nbs,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " block rowlength %" PetscInt_FMT,i,nnz[i],nbs); 1585 } 1586 } 1587 1588 B->ops->mult = MatMult_SeqSBAIJ_N; 1589 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1590 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1591 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1592 1593 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL)); 1594 if (!flg) { 1595 switch (bs) { 1596 case 1: 1597 B->ops->mult = MatMult_SeqSBAIJ_1; 1598 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1599 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1600 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1601 break; 1602 case 2: 1603 B->ops->mult = MatMult_SeqSBAIJ_2; 1604 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1605 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1606 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1607 break; 1608 case 3: 1609 B->ops->mult = MatMult_SeqSBAIJ_3; 1610 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1611 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1612 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1613 break; 1614 case 4: 1615 B->ops->mult = MatMult_SeqSBAIJ_4; 1616 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1617 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1618 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1619 break; 1620 case 5: 1621 B->ops->mult = MatMult_SeqSBAIJ_5; 1622 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1623 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1624 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1625 break; 1626 case 6: 1627 B->ops->mult = MatMult_SeqSBAIJ_6; 1628 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1629 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1630 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1631 break; 1632 case 7: 1633 B->ops->mult = MatMult_SeqSBAIJ_7; 1634 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1635 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1636 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1637 break; 1638 } 1639 } 1640 1641 b->mbs = mbs; 1642 b->nbs = nbs; 1643 if (!skipallocation) { 1644 if (!b->imax) { 1645 PetscCall(PetscMalloc2(mbs,&b->imax,mbs,&b->ilen)); 1646 1647 b->free_imax_ilen = PETSC_TRUE; 1648 1649 PetscCall(PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt))); 1650 } 1651 if (!nnz) { 1652 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1653 else if (nz <= 0) nz = 1; 1654 nz = PetscMin(nbs,nz); 1655 for (i=0; i<mbs; i++) b->imax[i] = nz; 1656 PetscCall(PetscIntMultError(nz,mbs,&nz)); 1657 } else { 1658 PetscInt64 nz64 = 0; 1659 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 1660 PetscCall(PetscIntCast(nz64,&nz)); 1661 } 1662 /* b->ilen will count nonzeros in each block row so far. */ 1663 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1664 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1665 1666 /* allocate the matrix space */ 1667 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 1668 PetscCall(PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i)); 1669 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)))); 1670 PetscCall(PetscArrayzero(b->a,nz*bs2)); 1671 PetscCall(PetscArrayzero(b->j,nz)); 1672 1673 b->singlemalloc = PETSC_TRUE; 1674 1675 /* pointer to beginning of each row */ 1676 b->i[0] = 0; 1677 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1678 1679 b->free_a = PETSC_TRUE; 1680 b->free_ij = PETSC_TRUE; 1681 } else { 1682 b->free_a = PETSC_FALSE; 1683 b->free_ij = PETSC_FALSE; 1684 } 1685 1686 b->bs2 = bs2; 1687 b->nz = 0; 1688 b->maxnz = nz; 1689 b->inew = NULL; 1690 b->jnew = NULL; 1691 b->anew = NULL; 1692 b->a2anew = NULL; 1693 b->permute = PETSC_FALSE; 1694 1695 B->was_assembled = PETSC_FALSE; 1696 B->assembled = PETSC_FALSE; 1697 if (realalloc) PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 1698 PetscFunctionReturn(0); 1699 } 1700 1701 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1702 { 1703 PetscInt i,j,m,nz,anz, nz_max=0,*nnz; 1704 PetscScalar *values=NULL; 1705 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1706 1707 PetscFunctionBegin; 1708 PetscCheck(bs >= 1,PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %" PetscInt_FMT,bs); 1709 PetscCall(PetscLayoutSetBlockSize(B->rmap,bs)); 1710 PetscCall(PetscLayoutSetBlockSize(B->cmap,bs)); 1711 PetscCall(PetscLayoutSetUp(B->rmap)); 1712 PetscCall(PetscLayoutSetUp(B->cmap)); 1713 PetscCall(PetscLayoutGetBlockSize(B->rmap,&bs)); 1714 m = B->rmap->n/bs; 1715 1716 PetscCheck(!ii[0],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %" PetscInt_FMT,ii[0]); 1717 PetscCall(PetscMalloc1(m+1,&nnz)); 1718 for (i=0; i<m; i++) { 1719 nz = ii[i+1] - ii[i]; 1720 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT,i,nz); 1721 anz = 0; 1722 for (j=0; j<nz; j++) { 1723 /* count only values on the diagonal or above */ 1724 if (jj[ii[i] + j] >= i) { 1725 anz = nz - j; 1726 break; 1727 } 1728 } 1729 nz_max = PetscMax(nz_max,anz); 1730 nnz[i] = anz; 1731 } 1732 PetscCall(MatSeqSBAIJSetPreallocation(B,bs,0,nnz)); 1733 PetscCall(PetscFree(nnz)); 1734 1735 values = (PetscScalar*)V; 1736 if (!values) { 1737 PetscCall(PetscCalloc1(bs*bs*nz_max,&values)); 1738 } 1739 for (i=0; i<m; i++) { 1740 PetscInt ncols = ii[i+1] - ii[i]; 1741 const PetscInt *icols = jj + ii[i]; 1742 if (!roworiented || bs == 1) { 1743 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1744 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES)); 1745 } else { 1746 for (j=0; j<ncols; j++) { 1747 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1748 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES)); 1749 } 1750 } 1751 } 1752 if (!V) PetscCall(PetscFree(values)); 1753 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1754 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1755 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 1756 PetscFunctionReturn(0); 1757 } 1758 1759 /* 1760 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1761 */ 1762 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1763 { 1764 PetscBool flg = PETSC_FALSE; 1765 PetscInt bs = B->rmap->bs; 1766 1767 PetscFunctionBegin; 1768 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL)); 1769 if (flg) bs = 8; 1770 1771 if (!natural) { 1772 switch (bs) { 1773 case 1: 1774 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1775 break; 1776 case 2: 1777 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1778 break; 1779 case 3: 1780 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1781 break; 1782 case 4: 1783 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1784 break; 1785 case 5: 1786 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1787 break; 1788 case 6: 1789 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1790 break; 1791 case 7: 1792 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1793 break; 1794 default: 1795 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1796 break; 1797 } 1798 } else { 1799 switch (bs) { 1800 case 1: 1801 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1802 break; 1803 case 2: 1804 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1805 break; 1806 case 3: 1807 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1808 break; 1809 case 4: 1810 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1811 break; 1812 case 5: 1813 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1814 break; 1815 case 6: 1816 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1817 break; 1818 case 7: 1819 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1820 break; 1821 default: 1822 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1823 break; 1824 } 1825 } 1826 PetscFunctionReturn(0); 1827 } 1828 1829 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*); 1830 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*); 1831 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type) 1832 { 1833 PetscFunctionBegin; 1834 *type = MATSOLVERPETSC; 1835 PetscFunctionReturn(0); 1836 } 1837 1838 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1839 { 1840 PetscInt n = A->rmap->n; 1841 1842 PetscFunctionBegin; 1843 #if defined(PETSC_USE_COMPLEX) 1844 PetscCheck(A->hermitian != PETSC_BOOL3_TRUE || A->symmetric == PETSC_BOOL3_TRUE || (ftype != MAT_FACTOR_CHOLESKY && ftype != MAT_FACTOR_ICC),PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian CHOLESKY or ICC Factor is not supported"); 1845 #endif 1846 1847 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 1848 PetscCall(MatSetSizes(*B,n,n,n,n)); 1849 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1850 PetscCall(MatSetType(*B,MATSEQSBAIJ)); 1851 PetscCall(MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL)); 1852 1853 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1854 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1855 PetscCall(PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY])); 1856 PetscCall(PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC])); 1857 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1858 1859 (*B)->factortype = ftype; 1860 (*B)->canuseordering = PETSC_TRUE; 1861 PetscCall(PetscFree((*B)->solvertype)); 1862 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype)); 1863 PetscCall(PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc)); 1864 PetscFunctionReturn(0); 1865 } 1866 1867 /*@C 1868 MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored 1869 1870 Not Collective 1871 1872 Input Parameter: 1873 . mat - a MATSEQSBAIJ matrix 1874 1875 Output Parameter: 1876 . array - pointer to the data 1877 1878 Level: intermediate 1879 1880 .seealso: `MatSeqSBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1881 @*/ 1882 PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array) 1883 { 1884 PetscFunctionBegin; 1885 PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array)); 1886 PetscFunctionReturn(0); 1887 } 1888 1889 /*@C 1890 MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray() 1891 1892 Not Collective 1893 1894 Input Parameters: 1895 + mat - a MATSEQSBAIJ matrix 1896 - array - pointer to the data 1897 1898 Level: intermediate 1899 1900 .seealso: `MatSeqSBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1901 @*/ 1902 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array) 1903 { 1904 PetscFunctionBegin; 1905 PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array)); 1906 PetscFunctionReturn(0); 1907 } 1908 1909 /*MC 1910 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1911 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1912 1913 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1914 can call MatSetOption(Mat, MAT_HERMITIAN). 1915 1916 Options Database Keys: 1917 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1918 1919 Notes: 1920 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1921 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1922 the options database -mat_ignore_lower_triangular false it will generate an error if you try to set a value in the lower triangular portion. 1923 1924 The number of rows in the matrix must be less than or equal to the number of columns 1925 1926 Level: beginner 1927 1928 .seealso: `MatCreateSeqSBAIJ()`, `MatType`, `MATMPISBAIJ` 1929 M*/ 1930 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1931 { 1932 Mat_SeqSBAIJ *b; 1933 PetscMPIInt size; 1934 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1935 1936 PetscFunctionBegin; 1937 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 1938 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1939 1940 PetscCall(PetscNewLog(B,&b)); 1941 B->data = (void*)b; 1942 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 1943 1944 B->ops->destroy = MatDestroy_SeqSBAIJ; 1945 B->ops->view = MatView_SeqSBAIJ; 1946 b->row = NULL; 1947 b->icol = NULL; 1948 b->reallocs = 0; 1949 b->saved_values = NULL; 1950 b->inode.limit = 5; 1951 b->inode.max_limit = 5; 1952 1953 b->roworiented = PETSC_TRUE; 1954 b->nonew = 0; 1955 b->diag = NULL; 1956 b->solve_work = NULL; 1957 b->mult_work = NULL; 1958 B->spptr = NULL; 1959 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1960 b->keepnonzeropattern = PETSC_FALSE; 1961 1962 b->inew = NULL; 1963 b->jnew = NULL; 1964 b->anew = NULL; 1965 b->a2anew = NULL; 1966 b->permute = PETSC_FALSE; 1967 1968 b->ignore_ltriangular = PETSC_TRUE; 1969 1970 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL)); 1971 1972 b->getrow_utriangular = PETSC_FALSE; 1973 1974 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL)); 1975 1976 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ)); 1977 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ)); 1978 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ)); 1979 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ)); 1980 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ)); 1981 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ)); 1982 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ)); 1983 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ)); 1984 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ)); 1985 #if defined(PETSC_HAVE_ELEMENTAL) 1986 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental)); 1987 #endif 1988 #if defined(PETSC_HAVE_SCALAPACK) 1989 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK)); 1990 #endif 1991 1992 B->symmetry_eternal = PETSC_TRUE; 1993 B->structural_symmetry_eternal = PETSC_TRUE; 1994 B->symmetric = PETSC_BOOL3_TRUE; 1995 B->structurally_symmetric = PETSC_BOOL3_TRUE; 1996 #if defined(PETSC_USE_COMPLEX) 1997 B->hermitian = PETSC_BOOL3_FALSE; 1998 #else 1999 B->hermitian = PETSC_BOOL3_TRUE; 2000 #endif 2001 2002 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ)); 2003 2004 PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat"); 2005 PetscCall(PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL)); 2006 if (no_unroll) { 2007 PetscCall(PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n")); 2008 } 2009 PetscCall(PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL)); 2010 if (no_inode) PetscCall(PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n")); 2011 PetscCall(PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL)); 2012 PetscOptionsEnd(); 2013 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 2014 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2015 PetscFunctionReturn(0); 2016 } 2017 2018 /*@C 2019 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2020 compressed row) format. For good matrix assembly performance the 2021 user should preallocate the matrix storage by setting the parameter nz 2022 (or the array nnz). By setting these parameters accurately, performance 2023 during matrix assembly can be increased by more than a factor of 50. 2024 2025 Collective on Mat 2026 2027 Input Parameters: 2028 + B - the symmetric matrix 2029 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2030 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2031 . nz - number of block nonzeros per block row (same for all rows) 2032 - nnz - array containing the number of block nonzeros in the upper triangular plus 2033 diagonal portion of each block (possibly different for each block row) or NULL 2034 2035 Options Database Keys: 2036 + -mat_no_unroll - uses code that does not unroll the loops in the 2037 block calculations (much slower) 2038 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2039 2040 Level: intermediate 2041 2042 Notes: 2043 Specify the preallocated storage with either nz or nnz (not both). 2044 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2045 allocation. See Users-Manual: ch_mat for details. 2046 2047 You can call MatGetInfo() to get information on how effective the preallocation was; 2048 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2049 You can also run with the option -info and look for messages with the string 2050 malloc in them to see if additional memory allocation was needed. 2051 2052 If the nnz parameter is given then the nz parameter is ignored 2053 2054 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2055 @*/ 2056 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2057 { 2058 PetscFunctionBegin; 2059 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2060 PetscValidType(B,1); 2061 PetscValidLogicalCollectiveInt(B,bs,2); 2062 PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz)); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 /*@C 2067 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values 2068 2069 Input Parameters: 2070 + B - the matrix 2071 . bs - size of block, the blocks are ALWAYS square. 2072 . i - the indices into j for the start of each local row (starts with zero) 2073 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2074 - v - optional values in the matrix 2075 2076 Level: advanced 2077 2078 Notes: 2079 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2080 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2081 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2082 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2083 block column and the second index is over columns within a block. 2084 2085 Any entries below the diagonal are ignored 2086 2087 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 2088 and usually the numerical values as well 2089 2090 .seealso: `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValuesBlocked()`, `MatSeqSBAIJSetPreallocation()`, `MATSEQSBAIJ` 2091 @*/ 2092 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2093 { 2094 PetscFunctionBegin; 2095 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2096 PetscValidType(B,1); 2097 PetscValidLogicalCollectiveInt(B,bs,2); 2098 PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v)); 2099 PetscFunctionReturn(0); 2100 } 2101 2102 /*@C 2103 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2104 compressed row) format. For good matrix assembly performance the 2105 user should preallocate the matrix storage by setting the parameter nz 2106 (or the array nnz). By setting these parameters accurately, performance 2107 during matrix assembly can be increased by more than a factor of 50. 2108 2109 Collective 2110 2111 Input Parameters: 2112 + comm - MPI communicator, set to PETSC_COMM_SELF 2113 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2114 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2115 . m - number of rows, or number of columns 2116 . nz - number of block nonzeros per block row (same for all rows) 2117 - nnz - array containing the number of block nonzeros in the upper triangular plus 2118 diagonal portion of each block (possibly different for each block row) or NULL 2119 2120 Output Parameter: 2121 . A - the symmetric matrix 2122 2123 Options Database Keys: 2124 + -mat_no_unroll - uses code that does not unroll the loops in the 2125 block calculations (much slower) 2126 - -mat_block_size - size of the blocks to use 2127 2128 Level: intermediate 2129 2130 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2131 MatXXXXSetPreallocation() paradigm instead of this routine directly. 2132 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2133 2134 Notes: 2135 The number of rows and columns must be divisible by blocksize. 2136 This matrix type does not support complex Hermitian operation. 2137 2138 Specify the preallocated storage with either nz or nnz (not both). 2139 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2140 allocation. See Users-Manual: ch_mat for details. 2141 2142 If the nnz parameter is given then the nz parameter is ignored 2143 2144 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2145 @*/ 2146 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2147 { 2148 PetscFunctionBegin; 2149 PetscCall(MatCreate(comm,A)); 2150 PetscCall(MatSetSizes(*A,m,n,m,n)); 2151 PetscCall(MatSetType(*A,MATSEQSBAIJ)); 2152 PetscCall(MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz)); 2153 PetscFunctionReturn(0); 2154 } 2155 2156 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2157 { 2158 Mat C; 2159 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2160 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2161 2162 PetscFunctionBegin; 2163 PetscCheck(a->i[mbs] == nz,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2164 2165 *B = NULL; 2166 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2167 PetscCall(MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n)); 2168 PetscCall(MatSetBlockSizesFromMats(C,A,A)); 2169 PetscCall(MatSetType(C,MATSEQSBAIJ)); 2170 c = (Mat_SeqSBAIJ*)C->data; 2171 2172 C->preallocated = PETSC_TRUE; 2173 C->factortype = A->factortype; 2174 c->row = NULL; 2175 c->icol = NULL; 2176 c->saved_values = NULL; 2177 c->keepnonzeropattern = a->keepnonzeropattern; 2178 C->assembled = PETSC_TRUE; 2179 2180 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 2181 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 2182 c->bs2 = a->bs2; 2183 c->mbs = a->mbs; 2184 c->nbs = a->nbs; 2185 2186 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2187 c->imax = a->imax; 2188 c->ilen = a->ilen; 2189 c->free_imax_ilen = PETSC_FALSE; 2190 } else { 2191 PetscCall(PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen)); 2192 PetscCall(PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt))); 2193 for (i=0; i<mbs; i++) { 2194 c->imax[i] = a->imax[i]; 2195 c->ilen[i] = a->ilen[i]; 2196 } 2197 c->free_imax_ilen = PETSC_TRUE; 2198 } 2199 2200 /* allocate the matrix space */ 2201 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2202 PetscCall(PetscMalloc1(bs2*nz,&c->a)); 2203 PetscCall(PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar))); 2204 c->i = a->i; 2205 c->j = a->j; 2206 c->singlemalloc = PETSC_FALSE; 2207 c->free_a = PETSC_TRUE; 2208 c->free_ij = PETSC_FALSE; 2209 c->parent = A; 2210 PetscCall(PetscObjectReference((PetscObject)A)); 2211 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 2212 PetscCall(MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 2213 } else { 2214 PetscCall(PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i)); 2215 PetscCall(PetscArraycpy(c->i,a->i,mbs+1)); 2216 PetscCall(PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)))); 2217 c->singlemalloc = PETSC_TRUE; 2218 c->free_a = PETSC_TRUE; 2219 c->free_ij = PETSC_TRUE; 2220 } 2221 if (mbs > 0) { 2222 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2223 PetscCall(PetscArraycpy(c->j,a->j,nz)); 2224 } 2225 if (cpvalues == MAT_COPY_VALUES) { 2226 PetscCall(PetscArraycpy(c->a,a->a,bs2*nz)); 2227 } else { 2228 PetscCall(PetscArrayzero(c->a,bs2*nz)); 2229 } 2230 if (a->jshort) { 2231 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2232 /* if the parent matrix is reassembled, this child matrix will never notice */ 2233 PetscCall(PetscMalloc1(nz,&c->jshort)); 2234 PetscCall(PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short))); 2235 PetscCall(PetscArraycpy(c->jshort,a->jshort,nz)); 2236 2237 c->free_jshort = PETSC_TRUE; 2238 } 2239 } 2240 2241 c->roworiented = a->roworiented; 2242 c->nonew = a->nonew; 2243 2244 if (a->diag) { 2245 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2246 c->diag = a->diag; 2247 c->free_diag = PETSC_FALSE; 2248 } else { 2249 PetscCall(PetscMalloc1(mbs,&c->diag)); 2250 PetscCall(PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt))); 2251 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2252 c->free_diag = PETSC_TRUE; 2253 } 2254 } 2255 c->nz = a->nz; 2256 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2257 c->solve_work = NULL; 2258 c->mult_work = NULL; 2259 2260 *B = C; 2261 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 2262 PetscFunctionReturn(0); 2263 } 2264 2265 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2266 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2267 2268 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer) 2269 { 2270 PetscBool isbinary; 2271 2272 PetscFunctionBegin; 2273 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 2274 PetscCheck(isbinary,PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)mat)->type_name); 2275 PetscCall(MatLoad_SeqSBAIJ_Binary(mat,viewer)); 2276 PetscFunctionReturn(0); 2277 } 2278 2279 /*@ 2280 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2281 (upper triangular entries in CSR format) provided by the user. 2282 2283 Collective 2284 2285 Input Parameters: 2286 + comm - must be an MPI communicator of size 1 2287 . bs - size of block 2288 . m - number of rows 2289 . n - number of columns 2290 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that row block row of the matrix 2291 . j - column indices 2292 - a - matrix values 2293 2294 Output Parameter: 2295 . mat - the matrix 2296 2297 Level: advanced 2298 2299 Notes: 2300 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2301 once the matrix is destroyed 2302 2303 You cannot set new nonzero locations into this matrix, that will generate an error. 2304 2305 The i and j indices are 0 based 2306 2307 When block size is greater than 1 the matrix values must be stored using the SBAIJ storage format (see the SBAIJ code to determine this). For block size of 1 2308 it is the regular CSR format excluding the lower triangular elements. 2309 2310 .seealso: `MatCreate()`, `MatCreateSBAIJ()`, `MatCreateSeqSBAIJ()` 2311 2312 @*/ 2313 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 2314 { 2315 PetscInt ii; 2316 Mat_SeqSBAIJ *sbaij; 2317 2318 PetscFunctionBegin; 2319 PetscCheck(bs == 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %" PetscInt_FMT " > 1 is not supported yet",bs); 2320 PetscCheck(m == 0 || i[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2321 2322 PetscCall(MatCreate(comm,mat)); 2323 PetscCall(MatSetSizes(*mat,m,n,m,n)); 2324 PetscCall(MatSetType(*mat,MATSEQSBAIJ)); 2325 PetscCall(MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL)); 2326 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2327 PetscCall(PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen)); 2328 PetscCall(PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt))); 2329 2330 sbaij->i = i; 2331 sbaij->j = j; 2332 sbaij->a = a; 2333 2334 sbaij->singlemalloc = PETSC_FALSE; 2335 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2336 sbaij->free_a = PETSC_FALSE; 2337 sbaij->free_ij = PETSC_FALSE; 2338 sbaij->free_imax_ilen = PETSC_TRUE; 2339 2340 for (ii=0; ii<m; ii++) { 2341 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2342 PetscCheck(i[ii+1] >= i[ii],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT,ii,i[ii+1] - i[ii]); 2343 } 2344 if (PetscDefined(USE_DEBUG)) { 2345 for (ii=0; ii<sbaij->i[m]; ii++) { 2346 PetscCheck(j[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 2347 PetscCheck(j[ii] < n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index too large at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 2348 } 2349 } 2350 2351 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 2352 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 2353 PetscFunctionReturn(0); 2354 } 2355 2356 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2357 { 2358 PetscFunctionBegin; 2359 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat)); 2360 PetscFunctionReturn(0); 2361 } 2362