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) { 1348 PetscCall((*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES)); 1349 } 1350 } 1351 PetscCall(MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY)); 1352 PetscFunctionReturn(0); 1353 } 1354 1355 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a) 1356 { 1357 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data; 1358 1359 PetscFunctionBegin; 1360 if (!Y->preallocated || !aij->nz) { 1361 PetscCall(MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL)); 1362 } 1363 PetscCall(MatShift_Basic(Y,a)); 1364 PetscFunctionReturn(0); 1365 } 1366 1367 /* -------------------------------------------------------------------*/ 1368 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1369 MatGetRow_SeqSBAIJ, 1370 MatRestoreRow_SeqSBAIJ, 1371 MatMult_SeqSBAIJ_N, 1372 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1373 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1374 MatMultAdd_SeqSBAIJ_N, 1375 NULL, 1376 NULL, 1377 NULL, 1378 /* 10*/ NULL, 1379 NULL, 1380 MatCholeskyFactor_SeqSBAIJ, 1381 MatSOR_SeqSBAIJ, 1382 MatTranspose_SeqSBAIJ, 1383 /* 15*/ MatGetInfo_SeqSBAIJ, 1384 MatEqual_SeqSBAIJ, 1385 MatGetDiagonal_SeqSBAIJ, 1386 MatDiagonalScale_SeqSBAIJ, 1387 MatNorm_SeqSBAIJ, 1388 /* 20*/ NULL, 1389 MatAssemblyEnd_SeqSBAIJ, 1390 MatSetOption_SeqSBAIJ, 1391 MatZeroEntries_SeqSBAIJ, 1392 /* 24*/ NULL, 1393 NULL, 1394 NULL, 1395 NULL, 1396 NULL, 1397 /* 29*/ MatSetUp_SeqSBAIJ, 1398 NULL, 1399 NULL, 1400 NULL, 1401 NULL, 1402 /* 34*/ MatDuplicate_SeqSBAIJ, 1403 NULL, 1404 NULL, 1405 NULL, 1406 MatICCFactor_SeqSBAIJ, 1407 /* 39*/ MatAXPY_SeqSBAIJ, 1408 MatCreateSubMatrices_SeqSBAIJ, 1409 MatIncreaseOverlap_SeqSBAIJ, 1410 MatGetValues_SeqSBAIJ, 1411 MatCopy_SeqSBAIJ, 1412 /* 44*/ NULL, 1413 MatScale_SeqSBAIJ, 1414 MatShift_SeqSBAIJ, 1415 NULL, 1416 MatZeroRowsColumns_SeqSBAIJ, 1417 /* 49*/ NULL, 1418 MatGetRowIJ_SeqSBAIJ, 1419 MatRestoreRowIJ_SeqSBAIJ, 1420 NULL, 1421 NULL, 1422 /* 54*/ NULL, 1423 NULL, 1424 NULL, 1425 MatPermute_SeqSBAIJ, 1426 MatSetValuesBlocked_SeqSBAIJ, 1427 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1428 NULL, 1429 NULL, 1430 NULL, 1431 NULL, 1432 /* 64*/ NULL, 1433 NULL, 1434 NULL, 1435 NULL, 1436 NULL, 1437 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1438 NULL, 1439 MatConvert_MPISBAIJ_Basic, 1440 NULL, 1441 NULL, 1442 /* 74*/ NULL, 1443 NULL, 1444 NULL, 1445 NULL, 1446 NULL, 1447 /* 79*/ NULL, 1448 NULL, 1449 NULL, 1450 MatGetInertia_SeqSBAIJ, 1451 MatLoad_SeqSBAIJ, 1452 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1453 MatIsHermitian_SeqSBAIJ, 1454 MatIsStructurallySymmetric_SeqSBAIJ, 1455 NULL, 1456 NULL, 1457 /* 89*/ NULL, 1458 NULL, 1459 NULL, 1460 NULL, 1461 NULL, 1462 /* 94*/ NULL, 1463 NULL, 1464 NULL, 1465 NULL, 1466 NULL, 1467 /* 99*/ NULL, 1468 NULL, 1469 NULL, 1470 MatConjugate_SeqSBAIJ, 1471 NULL, 1472 /*104*/ NULL, 1473 MatRealPart_SeqSBAIJ, 1474 MatImaginaryPart_SeqSBAIJ, 1475 MatGetRowUpperTriangular_SeqSBAIJ, 1476 MatRestoreRowUpperTriangular_SeqSBAIJ, 1477 /*109*/ NULL, 1478 NULL, 1479 NULL, 1480 NULL, 1481 MatMissingDiagonal_SeqSBAIJ, 1482 /*114*/ NULL, 1483 NULL, 1484 NULL, 1485 NULL, 1486 NULL, 1487 /*119*/ NULL, 1488 NULL, 1489 NULL, 1490 NULL, 1491 NULL, 1492 /*124*/ NULL, 1493 NULL, 1494 NULL, 1495 NULL, 1496 NULL, 1497 /*129*/ NULL, 1498 NULL, 1499 NULL, 1500 NULL, 1501 NULL, 1502 /*134*/ NULL, 1503 NULL, 1504 NULL, 1505 NULL, 1506 NULL, 1507 /*139*/ MatSetBlockSizes_Default, 1508 NULL, 1509 NULL, 1510 NULL, 1511 NULL, 1512 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ, 1513 NULL, 1514 NULL, 1515 NULL, 1516 NULL, 1517 NULL, 1518 /*150*/ NULL 1519 }; 1520 1521 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1522 { 1523 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1524 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1525 1526 PetscFunctionBegin; 1527 PetscCheck(aij->nonew == 1,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1528 1529 /* allocate space for values if not already there */ 1530 if (!aij->saved_values) { 1531 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 1532 } 1533 1534 /* copy values over */ 1535 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 1536 PetscFunctionReturn(0); 1537 } 1538 1539 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1540 { 1541 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1542 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1543 1544 PetscFunctionBegin; 1545 PetscCheck(aij->nonew == 1,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1546 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1547 1548 /* copy values over */ 1549 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 1550 PetscFunctionReturn(0); 1551 } 1552 1553 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1554 { 1555 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1556 PetscInt i,mbs,nbs,bs2; 1557 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1558 1559 PetscFunctionBegin; 1560 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1561 1562 PetscCall(MatSetBlockSize(B,PetscAbs(bs))); 1563 PetscCall(PetscLayoutSetUp(B->rmap)); 1564 PetscCall(PetscLayoutSetUp(B->cmap)); 1565 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); 1566 PetscCall(PetscLayoutGetBlockSize(B->rmap,&bs)); 1567 1568 B->preallocated = PETSC_TRUE; 1569 1570 mbs = B->rmap->N/bs; 1571 nbs = B->cmap->n/bs; 1572 bs2 = bs*bs; 1573 1574 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"); 1575 1576 if (nz == MAT_SKIP_ALLOCATION) { 1577 skipallocation = PETSC_TRUE; 1578 nz = 0; 1579 } 1580 1581 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1582 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 1583 if (nnz) { 1584 for (i=0; i<mbs; i++) { 1585 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]); 1586 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); 1587 } 1588 } 1589 1590 B->ops->mult = MatMult_SeqSBAIJ_N; 1591 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1592 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1593 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1594 1595 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL)); 1596 if (!flg) { 1597 switch (bs) { 1598 case 1: 1599 B->ops->mult = MatMult_SeqSBAIJ_1; 1600 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1601 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1602 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1603 break; 1604 case 2: 1605 B->ops->mult = MatMult_SeqSBAIJ_2; 1606 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1607 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1608 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1609 break; 1610 case 3: 1611 B->ops->mult = MatMult_SeqSBAIJ_3; 1612 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1613 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1614 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1615 break; 1616 case 4: 1617 B->ops->mult = MatMult_SeqSBAIJ_4; 1618 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1619 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1620 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1621 break; 1622 case 5: 1623 B->ops->mult = MatMult_SeqSBAIJ_5; 1624 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1625 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1626 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1627 break; 1628 case 6: 1629 B->ops->mult = MatMult_SeqSBAIJ_6; 1630 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1631 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1632 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1633 break; 1634 case 7: 1635 B->ops->mult = MatMult_SeqSBAIJ_7; 1636 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1637 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1638 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1639 break; 1640 } 1641 } 1642 1643 b->mbs = mbs; 1644 b->nbs = nbs; 1645 if (!skipallocation) { 1646 if (!b->imax) { 1647 PetscCall(PetscMalloc2(mbs,&b->imax,mbs,&b->ilen)); 1648 1649 b->free_imax_ilen = PETSC_TRUE; 1650 1651 PetscCall(PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt))); 1652 } 1653 if (!nnz) { 1654 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1655 else if (nz <= 0) nz = 1; 1656 nz = PetscMin(nbs,nz); 1657 for (i=0; i<mbs; i++) b->imax[i] = nz; 1658 PetscCall(PetscIntMultError(nz,mbs,&nz)); 1659 } else { 1660 PetscInt64 nz64 = 0; 1661 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 1662 PetscCall(PetscIntCast(nz64,&nz)); 1663 } 1664 /* b->ilen will count nonzeros in each block row so far. */ 1665 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1666 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1667 1668 /* allocate the matrix space */ 1669 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 1670 PetscCall(PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i)); 1671 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)))); 1672 PetscCall(PetscArrayzero(b->a,nz*bs2)); 1673 PetscCall(PetscArrayzero(b->j,nz)); 1674 1675 b->singlemalloc = PETSC_TRUE; 1676 1677 /* pointer to beginning of each row */ 1678 b->i[0] = 0; 1679 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1680 1681 b->free_a = PETSC_TRUE; 1682 b->free_ij = PETSC_TRUE; 1683 } else { 1684 b->free_a = PETSC_FALSE; 1685 b->free_ij = PETSC_FALSE; 1686 } 1687 1688 b->bs2 = bs2; 1689 b->nz = 0; 1690 b->maxnz = nz; 1691 b->inew = NULL; 1692 b->jnew = NULL; 1693 b->anew = NULL; 1694 b->a2anew = NULL; 1695 b->permute = PETSC_FALSE; 1696 1697 B->was_assembled = PETSC_FALSE; 1698 B->assembled = PETSC_FALSE; 1699 if (realalloc) PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 1700 PetscFunctionReturn(0); 1701 } 1702 1703 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1704 { 1705 PetscInt i,j,m,nz,anz, nz_max=0,*nnz; 1706 PetscScalar *values=NULL; 1707 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1708 1709 PetscFunctionBegin; 1710 PetscCheck(bs >= 1,PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %" PetscInt_FMT,bs); 1711 PetscCall(PetscLayoutSetBlockSize(B->rmap,bs)); 1712 PetscCall(PetscLayoutSetBlockSize(B->cmap,bs)); 1713 PetscCall(PetscLayoutSetUp(B->rmap)); 1714 PetscCall(PetscLayoutSetUp(B->cmap)); 1715 PetscCall(PetscLayoutGetBlockSize(B->rmap,&bs)); 1716 m = B->rmap->n/bs; 1717 1718 PetscCheck(!ii[0],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %" PetscInt_FMT,ii[0]); 1719 PetscCall(PetscMalloc1(m+1,&nnz)); 1720 for (i=0; i<m; i++) { 1721 nz = ii[i+1] - ii[i]; 1722 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT,i,nz); 1723 anz = 0; 1724 for (j=0; j<nz; j++) { 1725 /* count only values on the diagonal or above */ 1726 if (jj[ii[i] + j] >= i) { 1727 anz = nz - j; 1728 break; 1729 } 1730 } 1731 nz_max = PetscMax(nz_max,anz); 1732 nnz[i] = anz; 1733 } 1734 PetscCall(MatSeqSBAIJSetPreallocation(B,bs,0,nnz)); 1735 PetscCall(PetscFree(nnz)); 1736 1737 values = (PetscScalar*)V; 1738 if (!values) { 1739 PetscCall(PetscCalloc1(bs*bs*nz_max,&values)); 1740 } 1741 for (i=0; i<m; i++) { 1742 PetscInt ncols = ii[i+1] - ii[i]; 1743 const PetscInt *icols = jj + ii[i]; 1744 if (!roworiented || bs == 1) { 1745 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1746 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES)); 1747 } else { 1748 for (j=0; j<ncols; j++) { 1749 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1750 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES)); 1751 } 1752 } 1753 } 1754 if (!V) PetscCall(PetscFree(values)); 1755 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1756 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1757 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 1758 PetscFunctionReturn(0); 1759 } 1760 1761 /* 1762 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1763 */ 1764 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1765 { 1766 PetscBool flg = PETSC_FALSE; 1767 PetscInt bs = B->rmap->bs; 1768 1769 PetscFunctionBegin; 1770 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL)); 1771 if (flg) bs = 8; 1772 1773 if (!natural) { 1774 switch (bs) { 1775 case 1: 1776 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1777 break; 1778 case 2: 1779 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1780 break; 1781 case 3: 1782 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1783 break; 1784 case 4: 1785 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1786 break; 1787 case 5: 1788 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1789 break; 1790 case 6: 1791 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1792 break; 1793 case 7: 1794 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1795 break; 1796 default: 1797 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1798 break; 1799 } 1800 } else { 1801 switch (bs) { 1802 case 1: 1803 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1804 break; 1805 case 2: 1806 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1807 break; 1808 case 3: 1809 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1810 break; 1811 case 4: 1812 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1813 break; 1814 case 5: 1815 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1816 break; 1817 case 6: 1818 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1819 break; 1820 case 7: 1821 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1822 break; 1823 default: 1824 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1825 break; 1826 } 1827 } 1828 PetscFunctionReturn(0); 1829 } 1830 1831 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*); 1832 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*); 1833 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type) 1834 { 1835 PetscFunctionBegin; 1836 *type = MATSOLVERPETSC; 1837 PetscFunctionReturn(0); 1838 } 1839 1840 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1841 { 1842 PetscInt n = A->rmap->n; 1843 1844 PetscFunctionBegin; 1845 #if defined(PETSC_USE_COMPLEX) 1846 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"); 1847 #endif 1848 1849 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 1850 PetscCall(MatSetSizes(*B,n,n,n,n)); 1851 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1852 PetscCall(MatSetType(*B,MATSEQSBAIJ)); 1853 PetscCall(MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL)); 1854 1855 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1856 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1857 PetscCall(PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY])); 1858 PetscCall(PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC])); 1859 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1860 1861 (*B)->factortype = ftype; 1862 (*B)->canuseordering = PETSC_TRUE; 1863 PetscCall(PetscFree((*B)->solvertype)); 1864 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype)); 1865 PetscCall(PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc)); 1866 PetscFunctionReturn(0); 1867 } 1868 1869 /*@C 1870 MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored 1871 1872 Not Collective 1873 1874 Input Parameter: 1875 . mat - a MATSEQSBAIJ matrix 1876 1877 Output Parameter: 1878 . array - pointer to the data 1879 1880 Level: intermediate 1881 1882 .seealso: `MatSeqSBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1883 @*/ 1884 PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array) 1885 { 1886 PetscFunctionBegin; 1887 PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array)); 1888 PetscFunctionReturn(0); 1889 } 1890 1891 /*@C 1892 MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray() 1893 1894 Not Collective 1895 1896 Input Parameters: 1897 + mat - a MATSEQSBAIJ matrix 1898 - array - pointer to the data 1899 1900 Level: intermediate 1901 1902 .seealso: `MatSeqSBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1903 @*/ 1904 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array) 1905 { 1906 PetscFunctionBegin; 1907 PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array)); 1908 PetscFunctionReturn(0); 1909 } 1910 1911 /*MC 1912 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1913 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1914 1915 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1916 can call MatSetOption(Mat, MAT_HERMITIAN). 1917 1918 Options Database Keys: 1919 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1920 1921 Notes: 1922 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1923 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1924 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. 1925 1926 The number of rows in the matrix must be less than or equal to the number of columns 1927 1928 Level: beginner 1929 1930 .seealso: `MatCreateSeqSBAIJ()`, `MatType`, `MATMPISBAIJ` 1931 M*/ 1932 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1933 { 1934 Mat_SeqSBAIJ *b; 1935 PetscMPIInt size; 1936 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1937 1938 PetscFunctionBegin; 1939 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 1940 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1941 1942 PetscCall(PetscNewLog(B,&b)); 1943 B->data = (void*)b; 1944 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 1945 1946 B->ops->destroy = MatDestroy_SeqSBAIJ; 1947 B->ops->view = MatView_SeqSBAIJ; 1948 b->row = NULL; 1949 b->icol = NULL; 1950 b->reallocs = 0; 1951 b->saved_values = NULL; 1952 b->inode.limit = 5; 1953 b->inode.max_limit = 5; 1954 1955 b->roworiented = PETSC_TRUE; 1956 b->nonew = 0; 1957 b->diag = NULL; 1958 b->solve_work = NULL; 1959 b->mult_work = NULL; 1960 B->spptr = NULL; 1961 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1962 b->keepnonzeropattern = PETSC_FALSE; 1963 1964 b->inew = NULL; 1965 b->jnew = NULL; 1966 b->anew = NULL; 1967 b->a2anew = NULL; 1968 b->permute = PETSC_FALSE; 1969 1970 b->ignore_ltriangular = PETSC_TRUE; 1971 1972 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL)); 1973 1974 b->getrow_utriangular = PETSC_FALSE; 1975 1976 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL)); 1977 1978 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ)); 1979 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ)); 1980 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ)); 1981 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ)); 1982 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ)); 1983 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ)); 1984 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ)); 1985 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ)); 1986 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ)); 1987 #if defined(PETSC_HAVE_ELEMENTAL) 1988 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental)); 1989 #endif 1990 #if defined(PETSC_HAVE_SCALAPACK) 1991 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK)); 1992 #endif 1993 1994 B->symmetry_eternal = PETSC_TRUE; 1995 B->structural_symmetry_eternal = PETSC_TRUE; 1996 B->symmetric = PETSC_BOOL3_TRUE; 1997 B->structurally_symmetric = PETSC_BOOL3_TRUE; 1998 #if defined(PETSC_USE_COMPLEX) 1999 B->hermitian = PETSC_BOOL3_FALSE; 2000 #else 2001 B->hermitian = PETSC_BOOL3_TRUE; 2002 #endif 2003 2004 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ)); 2005 2006 PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat"); 2007 PetscCall(PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL)); 2008 if (no_unroll) { 2009 PetscCall(PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n")); 2010 } 2011 PetscCall(PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL)); 2012 if (no_inode) PetscCall(PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n")); 2013 PetscCall(PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL)); 2014 PetscOptionsEnd(); 2015 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 2016 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2017 PetscFunctionReturn(0); 2018 } 2019 2020 /*@C 2021 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2022 compressed row) format. For good matrix assembly performance the 2023 user should preallocate the matrix storage by setting the parameter nz 2024 (or the array nnz). By setting these parameters accurately, performance 2025 during matrix assembly can be increased by more than a factor of 50. 2026 2027 Collective on Mat 2028 2029 Input Parameters: 2030 + B - the symmetric matrix 2031 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2032 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2033 . nz - number of block nonzeros per block row (same for all rows) 2034 - nnz - array containing the number of block nonzeros in the upper triangular plus 2035 diagonal portion of each block (possibly different for each block row) or NULL 2036 2037 Options Database Keys: 2038 + -mat_no_unroll - uses code that does not unroll the loops in the 2039 block calculations (much slower) 2040 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2041 2042 Level: intermediate 2043 2044 Notes: 2045 Specify the preallocated storage with either nz or nnz (not both). 2046 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2047 allocation. See Users-Manual: ch_mat for details. 2048 2049 You can call MatGetInfo() to get information on how effective the preallocation was; 2050 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2051 You can also run with the option -info and look for messages with the string 2052 malloc in them to see if additional memory allocation was needed. 2053 2054 If the nnz parameter is given then the nz parameter is ignored 2055 2056 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2057 @*/ 2058 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2059 { 2060 PetscFunctionBegin; 2061 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2062 PetscValidType(B,1); 2063 PetscValidLogicalCollectiveInt(B,bs,2); 2064 PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz)); 2065 PetscFunctionReturn(0); 2066 } 2067 2068 /*@C 2069 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values 2070 2071 Input Parameters: 2072 + B - the matrix 2073 . bs - size of block, the blocks are ALWAYS square. 2074 . i - the indices into j for the start of each local row (starts with zero) 2075 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2076 - v - optional values in the matrix 2077 2078 Level: advanced 2079 2080 Notes: 2081 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2082 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2083 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2084 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2085 block column and the second index is over columns within a block. 2086 2087 Any entries below the diagonal are ignored 2088 2089 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 2090 and usually the numerical values as well 2091 2092 .seealso: `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValuesBlocked()`, `MatSeqSBAIJSetPreallocation()`, `MATSEQSBAIJ` 2093 @*/ 2094 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2095 { 2096 PetscFunctionBegin; 2097 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2098 PetscValidType(B,1); 2099 PetscValidLogicalCollectiveInt(B,bs,2); 2100 PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v)); 2101 PetscFunctionReturn(0); 2102 } 2103 2104 /*@C 2105 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2106 compressed row) format. For good matrix assembly performance the 2107 user should preallocate the matrix storage by setting the parameter nz 2108 (or the array nnz). By setting these parameters accurately, performance 2109 during matrix assembly can be increased by more than a factor of 50. 2110 2111 Collective 2112 2113 Input Parameters: 2114 + comm - MPI communicator, set to PETSC_COMM_SELF 2115 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2116 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2117 . m - number of rows, or number of columns 2118 . nz - number of block nonzeros per block row (same for all rows) 2119 - nnz - array containing the number of block nonzeros in the upper triangular plus 2120 diagonal portion of each block (possibly different for each block row) or NULL 2121 2122 Output Parameter: 2123 . A - the symmetric matrix 2124 2125 Options Database Keys: 2126 + -mat_no_unroll - uses code that does not unroll the loops in the 2127 block calculations (much slower) 2128 - -mat_block_size - size of the blocks to use 2129 2130 Level: intermediate 2131 2132 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2133 MatXXXXSetPreallocation() paradigm instead of this routine directly. 2134 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2135 2136 Notes: 2137 The number of rows and columns must be divisible by blocksize. 2138 This matrix type does not support complex Hermitian operation. 2139 2140 Specify the preallocated storage with either nz or nnz (not both). 2141 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2142 allocation. See Users-Manual: ch_mat for details. 2143 2144 If the nnz parameter is given then the nz parameter is ignored 2145 2146 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2147 @*/ 2148 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2149 { 2150 PetscFunctionBegin; 2151 PetscCall(MatCreate(comm,A)); 2152 PetscCall(MatSetSizes(*A,m,n,m,n)); 2153 PetscCall(MatSetType(*A,MATSEQSBAIJ)); 2154 PetscCall(MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz)); 2155 PetscFunctionReturn(0); 2156 } 2157 2158 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2159 { 2160 Mat C; 2161 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2162 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2163 2164 PetscFunctionBegin; 2165 PetscCheck(a->i[mbs] == nz,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2166 2167 *B = NULL; 2168 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2169 PetscCall(MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n)); 2170 PetscCall(MatSetBlockSizesFromMats(C,A,A)); 2171 PetscCall(MatSetType(C,MATSEQSBAIJ)); 2172 c = (Mat_SeqSBAIJ*)C->data; 2173 2174 C->preallocated = PETSC_TRUE; 2175 C->factortype = A->factortype; 2176 c->row = NULL; 2177 c->icol = NULL; 2178 c->saved_values = NULL; 2179 c->keepnonzeropattern = a->keepnonzeropattern; 2180 C->assembled = PETSC_TRUE; 2181 2182 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 2183 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 2184 c->bs2 = a->bs2; 2185 c->mbs = a->mbs; 2186 c->nbs = a->nbs; 2187 2188 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2189 c->imax = a->imax; 2190 c->ilen = a->ilen; 2191 c->free_imax_ilen = PETSC_FALSE; 2192 } else { 2193 PetscCall(PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen)); 2194 PetscCall(PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt))); 2195 for (i=0; i<mbs; i++) { 2196 c->imax[i] = a->imax[i]; 2197 c->ilen[i] = a->ilen[i]; 2198 } 2199 c->free_imax_ilen = PETSC_TRUE; 2200 } 2201 2202 /* allocate the matrix space */ 2203 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2204 PetscCall(PetscMalloc1(bs2*nz,&c->a)); 2205 PetscCall(PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar))); 2206 c->i = a->i; 2207 c->j = a->j; 2208 c->singlemalloc = PETSC_FALSE; 2209 c->free_a = PETSC_TRUE; 2210 c->free_ij = PETSC_FALSE; 2211 c->parent = A; 2212 PetscCall(PetscObjectReference((PetscObject)A)); 2213 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 2214 PetscCall(MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 2215 } else { 2216 PetscCall(PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i)); 2217 PetscCall(PetscArraycpy(c->i,a->i,mbs+1)); 2218 PetscCall(PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)))); 2219 c->singlemalloc = PETSC_TRUE; 2220 c->free_a = PETSC_TRUE; 2221 c->free_ij = PETSC_TRUE; 2222 } 2223 if (mbs > 0) { 2224 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2225 PetscCall(PetscArraycpy(c->j,a->j,nz)); 2226 } 2227 if (cpvalues == MAT_COPY_VALUES) { 2228 PetscCall(PetscArraycpy(c->a,a->a,bs2*nz)); 2229 } else { 2230 PetscCall(PetscArrayzero(c->a,bs2*nz)); 2231 } 2232 if (a->jshort) { 2233 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2234 /* if the parent matrix is reassembled, this child matrix will never notice */ 2235 PetscCall(PetscMalloc1(nz,&c->jshort)); 2236 PetscCall(PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short))); 2237 PetscCall(PetscArraycpy(c->jshort,a->jshort,nz)); 2238 2239 c->free_jshort = PETSC_TRUE; 2240 } 2241 } 2242 2243 c->roworiented = a->roworiented; 2244 c->nonew = a->nonew; 2245 2246 if (a->diag) { 2247 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2248 c->diag = a->diag; 2249 c->free_diag = PETSC_FALSE; 2250 } else { 2251 PetscCall(PetscMalloc1(mbs,&c->diag)); 2252 PetscCall(PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt))); 2253 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2254 c->free_diag = PETSC_TRUE; 2255 } 2256 } 2257 c->nz = a->nz; 2258 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2259 c->solve_work = NULL; 2260 c->mult_work = NULL; 2261 2262 *B = C; 2263 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 2264 PetscFunctionReturn(0); 2265 } 2266 2267 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2268 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2269 2270 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer) 2271 { 2272 PetscBool isbinary; 2273 2274 PetscFunctionBegin; 2275 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 2276 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); 2277 PetscCall(MatLoad_SeqSBAIJ_Binary(mat,viewer)); 2278 PetscFunctionReturn(0); 2279 } 2280 2281 /*@ 2282 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2283 (upper triangular entries in CSR format) provided by the user. 2284 2285 Collective 2286 2287 Input Parameters: 2288 + comm - must be an MPI communicator of size 1 2289 . bs - size of block 2290 . m - number of rows 2291 . n - number of columns 2292 . 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 2293 . j - column indices 2294 - a - matrix values 2295 2296 Output Parameter: 2297 . mat - the matrix 2298 2299 Level: advanced 2300 2301 Notes: 2302 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2303 once the matrix is destroyed 2304 2305 You cannot set new nonzero locations into this matrix, that will generate an error. 2306 2307 The i and j indices are 0 based 2308 2309 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 2310 it is the regular CSR format excluding the lower triangular elements. 2311 2312 .seealso: `MatCreate()`, `MatCreateSBAIJ()`, `MatCreateSeqSBAIJ()` 2313 2314 @*/ 2315 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 2316 { 2317 PetscInt ii; 2318 Mat_SeqSBAIJ *sbaij; 2319 2320 PetscFunctionBegin; 2321 PetscCheck(bs == 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %" PetscInt_FMT " > 1 is not supported yet",bs); 2322 PetscCheck(m == 0 || i[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2323 2324 PetscCall(MatCreate(comm,mat)); 2325 PetscCall(MatSetSizes(*mat,m,n,m,n)); 2326 PetscCall(MatSetType(*mat,MATSEQSBAIJ)); 2327 PetscCall(MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL)); 2328 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2329 PetscCall(PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen)); 2330 PetscCall(PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt))); 2331 2332 sbaij->i = i; 2333 sbaij->j = j; 2334 sbaij->a = a; 2335 2336 sbaij->singlemalloc = PETSC_FALSE; 2337 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2338 sbaij->free_a = PETSC_FALSE; 2339 sbaij->free_ij = PETSC_FALSE; 2340 sbaij->free_imax_ilen = PETSC_TRUE; 2341 2342 for (ii=0; ii<m; ii++) { 2343 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2344 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]); 2345 } 2346 if (PetscDefined(USE_DEBUG)) { 2347 for (ii=0; ii<sbaij->i[m]; ii++) { 2348 PetscCheck(j[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 2349 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]); 2350 } 2351 } 2352 2353 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 2354 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 2355 PetscFunctionReturn(0); 2356 } 2357 2358 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2359 { 2360 PetscFunctionBegin; 2361 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat)); 2362 PetscFunctionReturn(0); 2363 } 2364