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