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_INITIAL_MATRIX) { 335 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,B)); 336 } else if (reuse == MAT_REUSE_MATRIX) { 337 PetscCall(MatCopy(A,*B,SAME_NONZERO_PATTERN)); 338 } 339 PetscFunctionReturn(0); 340 } 341 342 PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 343 { 344 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 345 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 346 PetscViewerFormat format; 347 PetscInt *diag; 348 349 PetscFunctionBegin; 350 PetscCall(PetscViewerGetFormat(viewer,&format)); 351 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 352 PetscCall(PetscViewerASCIIPrintf(viewer," block size is %" PetscInt_FMT "\n",bs)); 353 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 354 Mat aij; 355 const char *matname; 356 357 if (A->factortype && bs>1) { 358 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n")); 359 PetscFunctionReturn(0); 360 } 361 PetscCall(MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij)); 362 PetscCall(PetscObjectGetName((PetscObject)A,&matname)); 363 PetscCall(PetscObjectSetName((PetscObject)aij,matname)); 364 PetscCall(MatView(aij,viewer)); 365 PetscCall(MatDestroy(&aij)); 366 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 367 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 368 for (i=0; i<a->mbs; i++) { 369 for (j=0; j<bs; j++) { 370 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i*bs+j)); 371 for (k=a->i[i]; k<a->i[i+1]; k++) { 372 for (l=0; l<bs; l++) { 373 #if defined(PETSC_USE_COMPLEX) 374 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 375 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k]+l, 376 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 377 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 378 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k]+l, 379 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 380 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 381 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]))); 382 } 383 #else 384 if (a->a[bs2*k + l*bs + j] != 0.0) { 385 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j])); 386 } 387 #endif 388 } 389 } 390 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 391 } 392 } 393 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 394 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 395 PetscFunctionReturn(0); 396 } else { 397 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 398 if (A->factortype) { /* for factored matrix */ 399 PetscCheck(bs<=1,PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet"); 400 401 diag=a->diag; 402 for (i=0; i<a->mbs; i++) { /* for row block i */ 403 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 404 /* diagonal entry */ 405 #if defined(PETSC_USE_COMPLEX) 406 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 407 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]]))); 408 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 409 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]]))); 410 } else { 411 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]))); 412 } 413 #else 414 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]))); 415 #endif 416 /* off-diagonal entries */ 417 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 418 #if defined(PETSC_USE_COMPLEX) 419 if (PetscImaginaryPart(a->a[k]) > 0.0) { 420 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]))); 421 } else if (PetscImaginaryPart(a->a[k]) < 0.0) { 422 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]))); 423 } else { 424 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]))); 425 } 426 #else 427 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[k],(double)a->a[k])); 428 #endif 429 } 430 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 431 } 432 433 } else { /* for non-factored matrix */ 434 for (i=0; i<a->mbs; i++) { /* for row block i */ 435 for (j=0; j<bs; j++) { /* for row bs*i + j */ 436 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i*bs+j)); 437 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */ 438 for (l=0; l<bs; l++) { /* for column */ 439 #if defined(PETSC_USE_COMPLEX) 440 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 441 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k]+l, 442 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 443 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 444 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k]+l, 445 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 446 } else { 447 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]))); 448 } 449 #else 450 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j])); 451 #endif 452 } 453 } 454 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 455 } 456 } 457 } 458 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 459 } 460 PetscCall(PetscViewerFlush(viewer)); 461 PetscFunctionReturn(0); 462 } 463 464 #include <petscdraw.h> 465 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 466 { 467 Mat A = (Mat) Aa; 468 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 469 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 470 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 471 MatScalar *aa; 472 PetscViewer viewer; 473 474 PetscFunctionBegin; 475 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 476 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 477 478 /* loop over matrix elements drawing boxes */ 479 480 PetscDrawCollectiveBegin(draw); 481 PetscCall(PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric")); 482 /* Blue for negative, Cyan for zero and Red for positive */ 483 color = PETSC_DRAW_BLUE; 484 for (i=0,row=0; i<mbs; i++,row+=bs) { 485 for (j=a->i[i]; j<a->i[i+1]; j++) { 486 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 487 x_l = a->j[j]*bs; x_r = x_l + 1.0; 488 aa = a->a + j*bs2; 489 for (k=0; k<bs; k++) { 490 for (l=0; l<bs; l++) { 491 if (PetscRealPart(*aa++) >= 0.) continue; 492 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 493 } 494 } 495 } 496 } 497 color = PETSC_DRAW_CYAN; 498 for (i=0,row=0; i<mbs; i++,row+=bs) { 499 for (j=a->i[i]; j<a->i[i+1]; j++) { 500 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 501 x_l = a->j[j]*bs; x_r = x_l + 1.0; 502 aa = a->a + j*bs2; 503 for (k=0; k<bs; k++) { 504 for (l=0; l<bs; l++) { 505 if (PetscRealPart(*aa++) != 0.) continue; 506 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 507 } 508 } 509 } 510 } 511 color = PETSC_DRAW_RED; 512 for (i=0,row=0; i<mbs; i++,row+=bs) { 513 for (j=a->i[i]; j<a->i[i+1]; j++) { 514 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 515 x_l = a->j[j]*bs; x_r = x_l + 1.0; 516 aa = a->a + j*bs2; 517 for (k=0; k<bs; k++) { 518 for (l=0; l<bs; l++) { 519 if (PetscRealPart(*aa++) <= 0.) continue; 520 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 521 } 522 } 523 } 524 } 525 PetscDrawCollectiveEnd(draw); 526 PetscFunctionReturn(0); 527 } 528 529 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 530 { 531 PetscReal xl,yl,xr,yr,w,h; 532 PetscDraw draw; 533 PetscBool isnull; 534 535 PetscFunctionBegin; 536 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 537 PetscCall(PetscDrawIsNull(draw,&isnull)); 538 if (isnull) PetscFunctionReturn(0); 539 540 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 541 xr += w; yr += h; xl = -w; yl = -h; 542 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 543 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 544 PetscCall(PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A)); 545 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 546 PetscCall(PetscDrawSave(draw)); 547 PetscFunctionReturn(0); 548 } 549 550 /* Used for both MPIBAIJ and MPISBAIJ matrices */ 551 #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary 552 553 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 554 { 555 PetscBool iascii,isbinary,isdraw; 556 557 PetscFunctionBegin; 558 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 559 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 560 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 561 if (iascii) { 562 PetscCall(MatView_SeqSBAIJ_ASCII(A,viewer)); 563 } else if (isbinary) { 564 PetscCall(MatView_SeqSBAIJ_Binary(A,viewer)); 565 } else if (isdraw) { 566 PetscCall(MatView_SeqSBAIJ_Draw(A,viewer)); 567 } else { 568 Mat B; 569 const char *matname; 570 PetscCall(MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B)); 571 PetscCall(PetscObjectGetName((PetscObject)A,&matname)); 572 PetscCall(PetscObjectSetName((PetscObject)B,matname)); 573 PetscCall(MatView(B,viewer)); 574 PetscCall(MatDestroy(&B)); 575 } 576 PetscFunctionReturn(0); 577 } 578 579 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 580 { 581 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 582 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 583 PetscInt *ai = a->i,*ailen = a->ilen; 584 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 585 MatScalar *ap,*aa = a->a; 586 587 PetscFunctionBegin; 588 for (k=0; k<m; k++) { /* loop over rows */ 589 row = im[k]; brow = row/bs; 590 if (row < 0) {v += n; continue;} /* negative row */ 591 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); 592 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 593 nrow = ailen[brow]; 594 for (l=0; l<n; l++) { /* loop over columns */ 595 if (in[l] < 0) {v++; continue;} /* negative column */ 596 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); 597 col = in[l]; 598 bcol = col/bs; 599 cidx = col%bs; 600 ridx = row%bs; 601 high = nrow; 602 low = 0; /* assume unsorted */ 603 while (high-low > 5) { 604 t = (low+high)/2; 605 if (rp[t] > bcol) high = t; 606 else low = t; 607 } 608 for (i=low; i<high; i++) { 609 if (rp[i] > bcol) break; 610 if (rp[i] == bcol) { 611 *v++ = ap[bs2*i+bs*cidx+ridx]; 612 goto finished; 613 } 614 } 615 *v++ = 0.0; 616 finished:; 617 } 618 } 619 PetscFunctionReturn(0); 620 } 621 622 PetscErrorCode MatPermute_SeqSBAIJ(Mat A,IS rowp,IS colp,Mat *B) 623 { 624 Mat C; 625 626 PetscFunctionBegin; 627 PetscCall(MatConvert(A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&C)); 628 PetscCall(MatPermute(C,rowp,colp,B)); 629 PetscCall(MatDestroy(&C)); 630 if (rowp == colp) { 631 PetscCall(MatConvert(*B,MATSEQSBAIJ,MAT_INPLACE_MATRIX,B)); 632 } 633 PetscFunctionReturn(0); 634 } 635 636 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 637 { 638 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 639 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 640 PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen; 641 PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 642 PetscBool roworiented=a->roworiented; 643 const PetscScalar *value = v; 644 MatScalar *ap,*aa = a->a,*bap; 645 646 PetscFunctionBegin; 647 if (roworiented) stepval = (n-1)*bs; 648 else stepval = (m-1)*bs; 649 650 for (k=0; k<m; k++) { /* loop over added rows */ 651 row = im[k]; 652 if (row < 0) continue; 653 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); 654 rp = aj + ai[row]; 655 ap = aa + bs2*ai[row]; 656 rmax = imax[row]; 657 nrow = ailen[row]; 658 low = 0; 659 high = nrow; 660 for (l=0; l<n; l++) { /* loop over added columns */ 661 if (in[l] < 0) continue; 662 col = in[l]; 663 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); 664 if (col < row) { 665 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */ 666 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)"); 667 } 668 if (roworiented) value = v + k*(stepval+bs)*bs + l*bs; 669 else value = v + l*(stepval+bs)*bs + k*bs; 670 671 if (col <= lastcol) low = 0; 672 else high = nrow; 673 674 lastcol = col; 675 while (high-low > 7) { 676 t = (low+high)/2; 677 if (rp[t] > col) high = t; 678 else low = t; 679 } 680 for (i=low; i<high; i++) { 681 if (rp[i] > col) break; 682 if (rp[i] == col) { 683 bap = ap + bs2*i; 684 if (roworiented) { 685 if (is == ADD_VALUES) { 686 for (ii=0; ii<bs; ii++,value+=stepval) { 687 for (jj=ii; jj<bs2; jj+=bs) { 688 bap[jj] += *value++; 689 } 690 } 691 } else { 692 for (ii=0; ii<bs; ii++,value+=stepval) { 693 for (jj=ii; jj<bs2; jj+=bs) { 694 bap[jj] = *value++; 695 } 696 } 697 } 698 } else { 699 if (is == ADD_VALUES) { 700 for (ii=0; ii<bs; ii++,value+=stepval) { 701 for (jj=0; jj<bs; jj++) { 702 *bap++ += *value++; 703 } 704 } 705 } else { 706 for (ii=0; ii<bs; ii++,value+=stepval) { 707 for (jj=0; jj<bs; jj++) { 708 *bap++ = *value++; 709 } 710 } 711 } 712 } 713 goto noinsert2; 714 } 715 } 716 if (nonew == 1) goto noinsert2; 717 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); 718 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 719 N = nrow++ - 1; high++; 720 /* shift up all the later entries in this row */ 721 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 722 PetscCall(PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1))); 723 PetscCall(PetscArrayzero(ap+bs2*i,bs2)); 724 rp[i] = col; 725 bap = ap + bs2*i; 726 if (roworiented) { 727 for (ii=0; ii<bs; ii++,value+=stepval) { 728 for (jj=ii; jj<bs2; jj+=bs) { 729 bap[jj] = *value++; 730 } 731 } 732 } else { 733 for (ii=0; ii<bs; ii++,value+=stepval) { 734 for (jj=0; jj<bs; jj++) { 735 *bap++ = *value++; 736 } 737 } 738 } 739 noinsert2:; 740 low = i; 741 } 742 ailen[row] = nrow; 743 } 744 PetscFunctionReturn(0); 745 } 746 747 /* 748 This is not yet used 749 */ 750 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 751 { 752 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 753 const PetscInt *ai = a->i, *aj = a->j,*cols; 754 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 755 PetscBool flag; 756 757 PetscFunctionBegin; 758 PetscCall(PetscMalloc1(m,&ns)); 759 while (i < m) { 760 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 761 /* Limits the number of elements in a node to 'a->inode.limit' */ 762 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 763 nzy = ai[j+1] - ai[j]; 764 if (nzy != (nzx - j + i)) break; 765 PetscCall(PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag)); 766 if (!flag) break; 767 } 768 ns[node_count++] = blk_size; 769 770 i = j; 771 } 772 if (!a->inode.size && m && node_count > .9*m) { 773 PetscCall(PetscFree(ns)); 774 PetscCall(PetscInfo(A,"Found %" PetscInt_FMT " nodes out of %" PetscInt_FMT " rows. Not using Inode routines\n",node_count,m)); 775 } else { 776 a->inode.node_count = node_count; 777 778 PetscCall(PetscMalloc1(node_count,&a->inode.size)); 779 PetscCall(PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt))); 780 PetscCall(PetscArraycpy(a->inode.size,ns,node_count)); 781 PetscCall(PetscFree(ns)); 782 PetscCall(PetscInfo(A,"Found %" PetscInt_FMT " nodes of %" PetscInt_FMT ". Limit used: %" PetscInt_FMT ". Using Inode routines\n",node_count,m,a->inode.limit)); 783 784 /* count collections of adjacent columns in each inode */ 785 row = 0; 786 cnt = 0; 787 for (i=0; i<node_count; i++) { 788 cols = aj + ai[row] + a->inode.size[i]; 789 nz = ai[row+1] - ai[row] - a->inode.size[i]; 790 for (j=1; j<nz; j++) { 791 if (cols[j] != cols[j-1]+1) cnt++; 792 } 793 cnt++; 794 row += a->inode.size[i]; 795 } 796 PetscCall(PetscMalloc1(2*cnt,&counts)); 797 cnt = 0; 798 row = 0; 799 for (i=0; i<node_count; i++) { 800 cols = aj + ai[row] + a->inode.size[i]; 801 counts[2*cnt] = cols[0]; 802 nz = ai[row+1] - ai[row] - a->inode.size[i]; 803 cnt2 = 1; 804 for (j=1; j<nz; j++) { 805 if (cols[j] != cols[j-1]+1) { 806 counts[2*(cnt++)+1] = cnt2; 807 counts[2*cnt] = cols[j]; 808 cnt2 = 1; 809 } else cnt2++; 810 } 811 counts[2*(cnt++)+1] = cnt2; 812 row += a->inode.size[i]; 813 } 814 PetscCall(PetscIntView(2*cnt,counts,NULL)); 815 } 816 PetscFunctionReturn(0); 817 } 818 819 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 820 { 821 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 822 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 823 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 824 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 825 MatScalar *aa = a->a,*ap; 826 827 PetscFunctionBegin; 828 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 829 830 if (m) rmax = ailen[0]; 831 for (i=1; i<mbs; i++) { 832 /* move each row back by the amount of empty slots (fshift) before it*/ 833 fshift += imax[i-1] - ailen[i-1]; 834 rmax = PetscMax(rmax,ailen[i]); 835 if (fshift) { 836 ip = aj + ai[i]; 837 ap = aa + bs2*ai[i]; 838 N = ailen[i]; 839 PetscCall(PetscArraymove(ip-fshift,ip,N)); 840 PetscCall(PetscArraymove(ap-bs2*fshift,ap,bs2*N)); 841 } 842 ai[i] = ai[i-1] + ailen[i-1]; 843 } 844 if (mbs) { 845 fshift += imax[mbs-1] - ailen[mbs-1]; 846 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 847 } 848 /* reset ilen and imax for each row */ 849 for (i=0; i<mbs; i++) { 850 ailen[i] = imax[i] = ai[i+1] - ai[i]; 851 } 852 a->nz = ai[mbs]; 853 854 /* diagonals may have moved, reset it */ 855 if (a->diag) PetscCall(PetscArraycpy(a->diag,ai,mbs)); 856 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); 857 858 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)); 859 PetscCall(PetscInfo(A,"Number of mallocs during MatSetValues is %" PetscInt_FMT "\n",a->reallocs)); 860 PetscCall(PetscInfo(A,"Most nonzeros blocks in any row is %" PetscInt_FMT "\n",rmax)); 861 862 A->info.mallocs += a->reallocs; 863 a->reallocs = 0; 864 A->info.nz_unneeded = (PetscReal)fshift*bs2; 865 a->idiagvalid = PETSC_FALSE; 866 a->rmax = rmax; 867 868 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 869 if (a->jshort && a->free_jshort) { 870 /* when matrix data structure is changed, previous jshort must be replaced */ 871 PetscCall(PetscFree(a->jshort)); 872 } 873 PetscCall(PetscMalloc1(a->i[A->rmap->n],&a->jshort)); 874 PetscCall(PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short))); 875 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 876 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 877 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 878 a->free_jshort = PETSC_TRUE; 879 } 880 PetscFunctionReturn(0); 881 } 882 883 /* 884 This function returns an array of flags which indicate the locations of contiguous 885 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 886 then the resulting sizes = [3,1,1,3,1] corresponding to sets [(0,1,2),(3),(5),(6,7,8),(9)] 887 Assume: sizes should be long enough to hold all the values. 888 */ 889 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 890 { 891 PetscInt i,j,k,row; 892 PetscBool flg; 893 894 PetscFunctionBegin; 895 for (i=0,j=0; i<n; j++) { 896 row = idx[i]; 897 if (row%bs!=0) { /* Not the beginning of a block */ 898 sizes[j] = 1; 899 i++; 900 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 901 sizes[j] = 1; /* Also makes sure at least 'bs' values exist for next else */ 902 i++; 903 } else { /* Beginning of the block, so check if the complete block exists */ 904 flg = PETSC_TRUE; 905 for (k=1; k<bs; k++) { 906 if (row+k != idx[i+k]) { /* break in the block */ 907 flg = PETSC_FALSE; 908 break; 909 } 910 } 911 if (flg) { /* No break in the bs */ 912 sizes[j] = bs; 913 i += bs; 914 } else { 915 sizes[j] = 1; 916 i++; 917 } 918 } 919 } 920 *bs_max = j; 921 PetscFunctionReturn(0); 922 } 923 924 /* Only add/insert a(i,j) with i<=j (blocks). 925 Any a(i,j) with i>j input by user is ingored. 926 */ 927 928 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 929 { 930 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 931 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 932 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 933 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 934 PetscInt ridx,cidx,bs2=a->bs2; 935 MatScalar *ap,value,*aa=a->a,*bap; 936 937 PetscFunctionBegin; 938 for (k=0; k<m; k++) { /* loop over added rows */ 939 row = im[k]; /* row number */ 940 brow = row/bs; /* block row number */ 941 if (row < 0) continue; 942 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); 943 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 944 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 945 rmax = imax[brow]; /* maximum space allocated for this row */ 946 nrow = ailen[brow]; /* actual length of this row */ 947 low = 0; 948 high = nrow; 949 for (l=0; l<n; l++) { /* loop over added columns */ 950 if (in[l] < 0) continue; 951 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); 952 col = in[l]; 953 bcol = col/bs; /* block col number */ 954 955 if (brow > bcol) { 956 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */ 957 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)"); 958 } 959 960 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 961 if ((brow==bcol && ridx<=cidx) || (brow<bcol)) { 962 /* element value a(k,l) */ 963 if (roworiented) value = v[l + k*n]; 964 else value = v[k + l*m]; 965 966 /* move pointer bap to a(k,l) quickly and add/insert value */ 967 if (col <= lastcol) low = 0; 968 else high = nrow; 969 970 lastcol = col; 971 while (high-low > 7) { 972 t = (low+high)/2; 973 if (rp[t] > bcol) high = t; 974 else low = t; 975 } 976 for (i=low; i<high; i++) { 977 if (rp[i] > bcol) break; 978 if (rp[i] == bcol) { 979 bap = ap + bs2*i + bs*cidx + ridx; 980 if (is == ADD_VALUES) *bap += value; 981 else *bap = value; 982 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 983 if (brow == bcol && ridx < cidx) { 984 bap = ap + bs2*i + bs*ridx + cidx; 985 if (is == ADD_VALUES) *bap += value; 986 else *bap = value; 987 } 988 goto noinsert1; 989 } 990 } 991 992 if (nonew == 1) goto noinsert1; 993 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col); 994 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 995 996 N = nrow++ - 1; high++; 997 /* shift up all the later entries in this row */ 998 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 999 PetscCall(PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1))); 1000 PetscCall(PetscArrayzero(ap+bs2*i,bs2)); 1001 rp[i] = bcol; 1002 ap[bs2*i + bs*cidx + ridx] = value; 1003 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1004 if (brow == bcol && ridx < cidx) { 1005 ap[bs2*i + bs*ridx + cidx] = value; 1006 } 1007 A->nonzerostate++; 1008 noinsert1:; 1009 low = i; 1010 } 1011 } /* end of loop over added columns */ 1012 ailen[brow] = nrow; 1013 } /* end of loop over added rows */ 1014 PetscFunctionReturn(0); 1015 } 1016 1017 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1018 { 1019 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1020 Mat outA; 1021 PetscBool row_identity; 1022 1023 PetscFunctionBegin; 1024 PetscCheck(info->levels == 0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1025 PetscCall(ISIdentity(row,&row_identity)); 1026 PetscCheck(row_identity,PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1027 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()! */ 1028 1029 outA = inA; 1030 inA->factortype = MAT_FACTOR_ICC; 1031 PetscCall(PetscFree(inA->solvertype)); 1032 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype)); 1033 1034 PetscCall(MatMarkDiagonal_SeqSBAIJ(inA)); 1035 PetscCall(MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity)); 1036 1037 PetscCall(PetscObjectReference((PetscObject)row)); 1038 PetscCall(ISDestroy(&a->row)); 1039 a->row = row; 1040 PetscCall(PetscObjectReference((PetscObject)row)); 1041 PetscCall(ISDestroy(&a->col)); 1042 a->col = row; 1043 1044 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1045 if (a->icol) PetscCall(ISInvertPermutation(row,PETSC_DECIDE, &a->icol)); 1046 PetscCall(PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol)); 1047 1048 if (!a->solve_work) { 1049 PetscCall(PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work)); 1050 PetscCall(PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar))); 1051 } 1052 1053 PetscCall(MatCholeskyFactorNumeric(outA,inA,info)); 1054 PetscFunctionReturn(0); 1055 } 1056 1057 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1058 { 1059 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data; 1060 PetscInt i,nz,n; 1061 1062 PetscFunctionBegin; 1063 nz = baij->maxnz; 1064 n = mat->cmap->n; 1065 for (i=0; i<nz; i++) baij->j[i] = indices[i]; 1066 1067 baij->nz = nz; 1068 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i]; 1069 1070 PetscCall(MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 1071 PetscFunctionReturn(0); 1072 } 1073 1074 /*@ 1075 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1076 in the matrix. 1077 1078 Input Parameters: 1079 + mat - the SeqSBAIJ matrix 1080 - indices - the column indices 1081 1082 Level: advanced 1083 1084 Notes: 1085 This can be called if you have precomputed the nonzero structure of the 1086 matrix and want to provide it to the matrix object to improve the performance 1087 of the MatSetValues() operation. 1088 1089 You MUST have set the correct numbers of nonzeros per row in the call to 1090 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1091 1092 MUST be called before any calls to MatSetValues() 1093 1094 .seealso: `MatCreateSeqSBAIJ` 1095 @*/ 1096 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1097 { 1098 PetscFunctionBegin; 1099 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1100 PetscValidIntPointer(indices,2); 1101 PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices)); 1102 PetscFunctionReturn(0); 1103 } 1104 1105 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1106 { 1107 PetscBool isbaij; 1108 1109 PetscFunctionBegin; 1110 PetscCall(PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"")); 1111 PetscCheck(isbaij,PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name); 1112 /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */ 1113 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 1114 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1115 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1116 1117 PetscCheck(a->i[a->mbs] == b->i[b->mbs],PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1118 PetscCheck(a->mbs == b->mbs,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different"); 1119 PetscCheck(a->bs2 == b->bs2,PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size"); 1120 PetscCall(PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs])); 1121 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 1122 } else { 1123 PetscCall(MatGetRowUpperTriangular(A)); 1124 PetscCall(MatCopy_Basic(A,B,str)); 1125 PetscCall(MatRestoreRowUpperTriangular(A)); 1126 } 1127 PetscFunctionReturn(0); 1128 } 1129 1130 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) 1131 { 1132 PetscFunctionBegin; 1133 PetscCall(MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL)); 1134 PetscFunctionReturn(0); 1135 } 1136 1137 static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1138 { 1139 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1140 1141 PetscFunctionBegin; 1142 *array = a->a; 1143 PetscFunctionReturn(0); 1144 } 1145 1146 static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1147 { 1148 PetscFunctionBegin; 1149 *array = NULL; 1150 PetscFunctionReturn(0); 1151 } 1152 1153 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz) 1154 { 1155 PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs; 1156 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data; 1157 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data; 1158 1159 PetscFunctionBegin; 1160 /* Set the number of nonzeros in the new matrix */ 1161 PetscCall(MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz)); 1162 PetscFunctionReturn(0); 1163 } 1164 1165 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1166 { 1167 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data; 1168 PetscInt bs=Y->rmap->bs,bs2=bs*bs; 1169 PetscBLASInt one = 1; 1170 1171 PetscFunctionBegin; 1172 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 1173 PetscBool e = x->nz == y->nz && x->mbs == y->mbs ? PETSC_TRUE : PETSC_FALSE; 1174 if (e) { 1175 PetscCall(PetscArraycmp(x->i,y->i,x->mbs+1,&e)); 1176 if (e) { 1177 PetscCall(PetscArraycmp(x->j,y->j,x->i[x->mbs],&e)); 1178 if (e) str = SAME_NONZERO_PATTERN; 1179 } 1180 } 1181 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatStructure is not SAME_NONZERO_PATTERN"); 1182 } 1183 if (str == SAME_NONZERO_PATTERN) { 1184 PetscScalar alpha = a; 1185 PetscBLASInt bnz; 1186 PetscCall(PetscBLASIntCast(x->nz*bs2,&bnz)); 1187 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 1188 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 1189 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1190 PetscCall(MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE)); 1191 PetscCall(MatAXPY_Basic(Y,a,X,str)); 1192 PetscCall(MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE)); 1193 } else { 1194 Mat B; 1195 PetscInt *nnz; 1196 PetscCheck(bs == X->rmap->bs,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size"); 1197 PetscCall(MatGetRowUpperTriangular(X)); 1198 PetscCall(MatGetRowUpperTriangular(Y)); 1199 PetscCall(PetscMalloc1(Y->rmap->N,&nnz)); 1200 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),&B)); 1201 PetscCall(PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name)); 1202 PetscCall(MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N)); 1203 PetscCall(MatSetBlockSizesFromMats(B,Y,Y)); 1204 PetscCall(MatSetType(B,((PetscObject)Y)->type_name)); 1205 PetscCall(MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz)); 1206 PetscCall(MatSeqSBAIJSetPreallocation(B,bs,0,nnz)); 1207 1208 PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 1209 1210 PetscCall(MatHeaderMerge(Y,&B)); 1211 PetscCall(PetscFree(nnz)); 1212 PetscCall(MatRestoreRowUpperTriangular(X)); 1213 PetscCall(MatRestoreRowUpperTriangular(Y)); 1214 } 1215 PetscFunctionReturn(0); 1216 } 1217 1218 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1219 { 1220 PetscFunctionBegin; 1221 *flg = PETSC_TRUE; 1222 PetscFunctionReturn(0); 1223 } 1224 1225 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1226 { 1227 PetscFunctionBegin; 1228 *flg = PETSC_TRUE; 1229 PetscFunctionReturn(0); 1230 } 1231 1232 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1233 { 1234 PetscFunctionBegin; 1235 *flg = PETSC_FALSE; 1236 PetscFunctionReturn(0); 1237 } 1238 1239 PetscErrorCode MatConjugate_SeqSBAIJ(Mat A) 1240 { 1241 #if defined(PETSC_USE_COMPLEX) 1242 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1243 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1244 MatScalar *aa = a->a; 1245 1246 PetscFunctionBegin; 1247 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1248 #else 1249 PetscFunctionBegin; 1250 #endif 1251 PetscFunctionReturn(0); 1252 } 1253 1254 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1255 { 1256 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1257 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1258 MatScalar *aa = a->a; 1259 1260 PetscFunctionBegin; 1261 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1262 PetscFunctionReturn(0); 1263 } 1264 1265 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1266 { 1267 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1268 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1269 MatScalar *aa = a->a; 1270 1271 PetscFunctionBegin; 1272 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1273 PetscFunctionReturn(0); 1274 } 1275 1276 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1277 { 1278 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1279 PetscInt i,j,k,count; 1280 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col; 1281 PetscScalar zero = 0.0; 1282 MatScalar *aa; 1283 const PetscScalar *xx; 1284 PetscScalar *bb; 1285 PetscBool *zeroed,vecs = PETSC_FALSE; 1286 1287 PetscFunctionBegin; 1288 /* fix right hand side if needed */ 1289 if (x && b) { 1290 PetscCall(VecGetArrayRead(x,&xx)); 1291 PetscCall(VecGetArray(b,&bb)); 1292 vecs = PETSC_TRUE; 1293 } 1294 1295 /* zero the columns */ 1296 PetscCall(PetscCalloc1(A->rmap->n,&zeroed)); 1297 for (i=0; i<is_n; i++) { 1298 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]); 1299 zeroed[is_idx[i]] = PETSC_TRUE; 1300 } 1301 if (vecs) { 1302 for (i=0; i<A->rmap->N; i++) { 1303 row = i/bs; 1304 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1305 for (k=0; k<bs; k++) { 1306 col = bs*baij->j[j] + k; 1307 if (col <= i) continue; 1308 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1309 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1310 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1311 } 1312 } 1313 } 1314 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1315 } 1316 1317 for (i=0; i<A->rmap->N; i++) { 1318 if (!zeroed[i]) { 1319 row = i/bs; 1320 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1321 for (k=0; k<bs; k++) { 1322 col = bs*baij->j[j] + k; 1323 if (zeroed[col]) { 1324 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1325 aa[0] = 0.0; 1326 } 1327 } 1328 } 1329 } 1330 } 1331 PetscCall(PetscFree(zeroed)); 1332 if (vecs) { 1333 PetscCall(VecRestoreArrayRead(x,&xx)); 1334 PetscCall(VecRestoreArray(b,&bb)); 1335 } 1336 1337 /* zero the rows */ 1338 for (i=0; i<is_n; i++) { 1339 row = is_idx[i]; 1340 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1341 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1342 for (k=0; k<count; k++) { 1343 aa[0] = zero; 1344 aa += bs; 1345 } 1346 if (diag != 0.0) { 1347 PetscCall((*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES)); 1348 } 1349 } 1350 PetscCall(MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY)); 1351 PetscFunctionReturn(0); 1352 } 1353 1354 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a) 1355 { 1356 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data; 1357 1358 PetscFunctionBegin; 1359 if (!Y->preallocated || !aij->nz) { 1360 PetscCall(MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL)); 1361 } 1362 PetscCall(MatShift_Basic(Y,a)); 1363 PetscFunctionReturn(0); 1364 } 1365 1366 /* -------------------------------------------------------------------*/ 1367 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1368 MatGetRow_SeqSBAIJ, 1369 MatRestoreRow_SeqSBAIJ, 1370 MatMult_SeqSBAIJ_N, 1371 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1372 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1373 MatMultAdd_SeqSBAIJ_N, 1374 NULL, 1375 NULL, 1376 NULL, 1377 /* 10*/ NULL, 1378 NULL, 1379 MatCholeskyFactor_SeqSBAIJ, 1380 MatSOR_SeqSBAIJ, 1381 MatTranspose_SeqSBAIJ, 1382 /* 15*/ MatGetInfo_SeqSBAIJ, 1383 MatEqual_SeqSBAIJ, 1384 MatGetDiagonal_SeqSBAIJ, 1385 MatDiagonalScale_SeqSBAIJ, 1386 MatNorm_SeqSBAIJ, 1387 /* 20*/ NULL, 1388 MatAssemblyEnd_SeqSBAIJ, 1389 MatSetOption_SeqSBAIJ, 1390 MatZeroEntries_SeqSBAIJ, 1391 /* 24*/ NULL, 1392 NULL, 1393 NULL, 1394 NULL, 1395 NULL, 1396 /* 29*/ MatSetUp_SeqSBAIJ, 1397 NULL, 1398 NULL, 1399 NULL, 1400 NULL, 1401 /* 34*/ MatDuplicate_SeqSBAIJ, 1402 NULL, 1403 NULL, 1404 NULL, 1405 MatICCFactor_SeqSBAIJ, 1406 /* 39*/ MatAXPY_SeqSBAIJ, 1407 MatCreateSubMatrices_SeqSBAIJ, 1408 MatIncreaseOverlap_SeqSBAIJ, 1409 MatGetValues_SeqSBAIJ, 1410 MatCopy_SeqSBAIJ, 1411 /* 44*/ NULL, 1412 MatScale_SeqSBAIJ, 1413 MatShift_SeqSBAIJ, 1414 NULL, 1415 MatZeroRowsColumns_SeqSBAIJ, 1416 /* 49*/ NULL, 1417 MatGetRowIJ_SeqSBAIJ, 1418 MatRestoreRowIJ_SeqSBAIJ, 1419 NULL, 1420 NULL, 1421 /* 54*/ NULL, 1422 NULL, 1423 NULL, 1424 MatPermute_SeqSBAIJ, 1425 MatSetValuesBlocked_SeqSBAIJ, 1426 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1427 NULL, 1428 NULL, 1429 NULL, 1430 NULL, 1431 /* 64*/ NULL, 1432 NULL, 1433 NULL, 1434 NULL, 1435 NULL, 1436 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1437 NULL, 1438 MatConvert_MPISBAIJ_Basic, 1439 NULL, 1440 NULL, 1441 /* 74*/ NULL, 1442 NULL, 1443 NULL, 1444 NULL, 1445 NULL, 1446 /* 79*/ NULL, 1447 NULL, 1448 NULL, 1449 MatGetInertia_SeqSBAIJ, 1450 MatLoad_SeqSBAIJ, 1451 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1452 MatIsHermitian_SeqSBAIJ, 1453 MatIsStructurallySymmetric_SeqSBAIJ, 1454 NULL, 1455 NULL, 1456 /* 89*/ NULL, 1457 NULL, 1458 NULL, 1459 NULL, 1460 NULL, 1461 /* 94*/ NULL, 1462 NULL, 1463 NULL, 1464 NULL, 1465 NULL, 1466 /* 99*/ NULL, 1467 NULL, 1468 NULL, 1469 MatConjugate_SeqSBAIJ, 1470 NULL, 1471 /*104*/ NULL, 1472 MatRealPart_SeqSBAIJ, 1473 MatImaginaryPart_SeqSBAIJ, 1474 MatGetRowUpperTriangular_SeqSBAIJ, 1475 MatRestoreRowUpperTriangular_SeqSBAIJ, 1476 /*109*/ NULL, 1477 NULL, 1478 NULL, 1479 NULL, 1480 MatMissingDiagonal_SeqSBAIJ, 1481 /*114*/ NULL, 1482 NULL, 1483 NULL, 1484 NULL, 1485 NULL, 1486 /*119*/ NULL, 1487 NULL, 1488 NULL, 1489 NULL, 1490 NULL, 1491 /*124*/ NULL, 1492 NULL, 1493 NULL, 1494 NULL, 1495 NULL, 1496 /*129*/ NULL, 1497 NULL, 1498 NULL, 1499 NULL, 1500 NULL, 1501 /*134*/ NULL, 1502 NULL, 1503 NULL, 1504 NULL, 1505 NULL, 1506 /*139*/ MatSetBlockSizes_Default, 1507 NULL, 1508 NULL, 1509 NULL, 1510 NULL, 1511 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ, 1512 NULL, 1513 NULL, 1514 NULL, 1515 NULL, 1516 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