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_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_STRUCTURE_ONLY: 271 /* These options are handled directly by MatSetOption() */ 272 break; 273 case MAT_IGNORE_LOWER_TRIANGULAR: 274 a->ignore_ltriangular = flg; 275 break; 276 case MAT_ERROR_LOWER_TRIANGULAR: 277 a->ignore_ltriangular = flg; 278 break; 279 case MAT_GETROW_UPPERTRIANGULAR: 280 a->getrow_utriangular = flg; 281 break; 282 case MAT_SUBMAT_SINGLEIS: 283 break; 284 default: 285 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 286 } 287 PetscFunctionReturn(0); 288 } 289 290 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 291 { 292 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 293 294 PetscFunctionBegin; 295 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()"); 296 297 /* Get the upper triangular part of the row */ 298 PetscCall(MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a)); 299 PetscFunctionReturn(0); 300 } 301 302 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 303 { 304 PetscFunctionBegin; 305 if (nz) *nz = 0; 306 if (idx) PetscCall(PetscFree(*idx)); 307 if (v) PetscCall(PetscFree(*v)); 308 PetscFunctionReturn(0); 309 } 310 311 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 312 { 313 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 314 315 PetscFunctionBegin; 316 a->getrow_utriangular = PETSC_TRUE; 317 PetscFunctionReturn(0); 318 } 319 320 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 321 { 322 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 323 324 PetscFunctionBegin; 325 a->getrow_utriangular = PETSC_FALSE; 326 PetscFunctionReturn(0); 327 } 328 329 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B) 330 { 331 PetscFunctionBegin; 332 if (reuse == MAT_INITIAL_MATRIX) { 333 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,B)); 334 } else if (reuse == MAT_REUSE_MATRIX) { 335 PetscCall(MatCopy(A,*B,SAME_NONZERO_PATTERN)); 336 } 337 PetscFunctionReturn(0); 338 } 339 340 PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 341 { 342 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 343 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 344 PetscViewerFormat format; 345 PetscInt *diag; 346 347 PetscFunctionBegin; 348 PetscCall(PetscViewerGetFormat(viewer,&format)); 349 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 350 PetscCall(PetscViewerASCIIPrintf(viewer," block size is %" PetscInt_FMT "\n",bs)); 351 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 352 Mat aij; 353 const char *matname; 354 355 if (A->factortype && bs>1) { 356 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n")); 357 PetscFunctionReturn(0); 358 } 359 PetscCall(MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij)); 360 PetscCall(PetscObjectGetName((PetscObject)A,&matname)); 361 PetscCall(PetscObjectSetName((PetscObject)aij,matname)); 362 PetscCall(MatView(aij,viewer)); 363 PetscCall(MatDestroy(&aij)); 364 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 365 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 366 for (i=0; i<a->mbs; i++) { 367 for (j=0; j<bs; j++) { 368 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i*bs+j)); 369 for (k=a->i[i]; k<a->i[i+1]; k++) { 370 for (l=0; l<bs; l++) { 371 #if defined(PETSC_USE_COMPLEX) 372 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 (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 376 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k]+l, 377 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 378 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 379 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]))); 380 } 381 #else 382 if (a->a[bs2*k + l*bs + j] != 0.0) { 383 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j])); 384 } 385 #endif 386 } 387 } 388 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 389 } 390 } 391 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 392 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 393 PetscFunctionReturn(0); 394 } else { 395 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 396 if (A->factortype) { /* for factored matrix */ 397 PetscCheck(bs<=1,PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet"); 398 399 diag=a->diag; 400 for (i=0; i<a->mbs; i++) { /* for row block i */ 401 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 402 /* diagonal entry */ 403 #if defined(PETSC_USE_COMPLEX) 404 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 405 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]]))); 406 } else 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 { 409 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]))); 410 } 411 #else 412 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]))); 413 #endif 414 /* off-diagonal entries */ 415 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 416 #if defined(PETSC_USE_COMPLEX) 417 if (PetscImaginaryPart(a->a[k]) > 0.0) { 418 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]))); 419 } else 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 { 422 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]))); 423 } 424 #else 425 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[k],(double)a->a[k])); 426 #endif 427 } 428 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 429 } 430 431 } else { /* for non-factored matrix */ 432 for (i=0; i<a->mbs; i++) { /* for row block i */ 433 for (j=0; j<bs; j++) { /* for row bs*i + j */ 434 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i*bs+j)); 435 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */ 436 for (l=0; l<bs; l++) { /* for column */ 437 #if defined(PETSC_USE_COMPLEX) 438 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 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 442 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i) ",bs*a->j[k]+l, 443 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]))); 444 } else { 445 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]))); 446 } 447 #else 448 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j])); 449 #endif 450 } 451 } 452 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 453 } 454 } 455 } 456 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 457 } 458 PetscCall(PetscViewerFlush(viewer)); 459 PetscFunctionReturn(0); 460 } 461 462 #include <petscdraw.h> 463 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 464 { 465 Mat A = (Mat) Aa; 466 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 467 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 468 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 469 MatScalar *aa; 470 PetscViewer viewer; 471 472 PetscFunctionBegin; 473 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 474 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 475 476 /* loop over matrix elements drawing boxes */ 477 478 PetscDrawCollectiveBegin(draw); 479 PetscCall(PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric")); 480 /* Blue for negative, Cyan for zero and Red for positive */ 481 color = PETSC_DRAW_BLUE; 482 for (i=0,row=0; i<mbs; i++,row+=bs) { 483 for (j=a->i[i]; j<a->i[i+1]; j++) { 484 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 485 x_l = a->j[j]*bs; x_r = x_l + 1.0; 486 aa = a->a + j*bs2; 487 for (k=0; k<bs; k++) { 488 for (l=0; l<bs; l++) { 489 if (PetscRealPart(*aa++) >= 0.) continue; 490 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 491 } 492 } 493 } 494 } 495 color = PETSC_DRAW_CYAN; 496 for (i=0,row=0; i<mbs; i++,row+=bs) { 497 for (j=a->i[i]; j<a->i[i+1]; j++) { 498 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 499 x_l = a->j[j]*bs; x_r = x_l + 1.0; 500 aa = a->a + j*bs2; 501 for (k=0; k<bs; k++) { 502 for (l=0; l<bs; l++) { 503 if (PetscRealPart(*aa++) != 0.) continue; 504 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 505 } 506 } 507 } 508 } 509 color = PETSC_DRAW_RED; 510 for (i=0,row=0; i<mbs; i++,row+=bs) { 511 for (j=a->i[i]; j<a->i[i+1]; j++) { 512 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 513 x_l = a->j[j]*bs; x_r = x_l + 1.0; 514 aa = a->a + j*bs2; 515 for (k=0; k<bs; k++) { 516 for (l=0; l<bs; l++) { 517 if (PetscRealPart(*aa++) <= 0.) continue; 518 PetscCall(PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color)); 519 } 520 } 521 } 522 } 523 PetscDrawCollectiveEnd(draw); 524 PetscFunctionReturn(0); 525 } 526 527 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 528 { 529 PetscReal xl,yl,xr,yr,w,h; 530 PetscDraw draw; 531 PetscBool isnull; 532 533 PetscFunctionBegin; 534 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 535 PetscCall(PetscDrawIsNull(draw,&isnull)); 536 if (isnull) PetscFunctionReturn(0); 537 538 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 539 xr += w; yr += h; xl = -w; yl = -h; 540 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 541 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 542 PetscCall(PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A)); 543 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 544 PetscCall(PetscDrawSave(draw)); 545 PetscFunctionReturn(0); 546 } 547 548 /* Used for both MPIBAIJ and MPISBAIJ matrices */ 549 #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary 550 551 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 552 { 553 PetscBool iascii,isbinary,isdraw; 554 555 PetscFunctionBegin; 556 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 557 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 558 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 559 if (iascii) { 560 PetscCall(MatView_SeqSBAIJ_ASCII(A,viewer)); 561 } else if (isbinary) { 562 PetscCall(MatView_SeqSBAIJ_Binary(A,viewer)); 563 } else if (isdraw) { 564 PetscCall(MatView_SeqSBAIJ_Draw(A,viewer)); 565 } else { 566 Mat B; 567 const char *matname; 568 PetscCall(MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B)); 569 PetscCall(PetscObjectGetName((PetscObject)A,&matname)); 570 PetscCall(PetscObjectSetName((PetscObject)B,matname)); 571 PetscCall(MatView(B,viewer)); 572 PetscCall(MatDestroy(&B)); 573 } 574 PetscFunctionReturn(0); 575 } 576 577 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 578 { 579 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 580 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 581 PetscInt *ai = a->i,*ailen = a->ilen; 582 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 583 MatScalar *ap,*aa = a->a; 584 585 PetscFunctionBegin; 586 for (k=0; k<m; k++) { /* loop over rows */ 587 row = im[k]; brow = row/bs; 588 if (row < 0) {v += n; continue;} /* negative row */ 589 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); 590 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 591 nrow = ailen[brow]; 592 for (l=0; l<n; l++) { /* loop over columns */ 593 if (in[l] < 0) {v++; continue;} /* negative column */ 594 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); 595 col = in[l]; 596 bcol = col/bs; 597 cidx = col%bs; 598 ridx = row%bs; 599 high = nrow; 600 low = 0; /* assume unsorted */ 601 while (high-low > 5) { 602 t = (low+high)/2; 603 if (rp[t] > bcol) high = t; 604 else low = t; 605 } 606 for (i=low; i<high; i++) { 607 if (rp[i] > bcol) break; 608 if (rp[i] == bcol) { 609 *v++ = ap[bs2*i+bs*cidx+ridx]; 610 goto finished; 611 } 612 } 613 *v++ = 0.0; 614 finished:; 615 } 616 } 617 PetscFunctionReturn(0); 618 } 619 620 PetscErrorCode MatPermute_SeqSBAIJ(Mat A,IS rowp,IS colp,Mat *B) 621 { 622 Mat C; 623 624 PetscFunctionBegin; 625 PetscCall(MatConvert(A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&C)); 626 PetscCall(MatPermute(C,rowp,colp,B)); 627 PetscCall(MatDestroy(&C)); 628 if (rowp == colp) { 629 PetscCall(MatConvert(*B,MATSEQSBAIJ,MAT_INPLACE_MATRIX,B)); 630 } 631 PetscFunctionReturn(0); 632 } 633 634 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 635 { 636 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 637 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 638 PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen; 639 PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 640 PetscBool roworiented=a->roworiented; 641 const PetscScalar *value = v; 642 MatScalar *ap,*aa = a->a,*bap; 643 644 PetscFunctionBegin; 645 if (roworiented) stepval = (n-1)*bs; 646 else stepval = (m-1)*bs; 647 648 for (k=0; k<m; k++) { /* loop over added rows */ 649 row = im[k]; 650 if (row < 0) continue; 651 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); 652 rp = aj + ai[row]; 653 ap = aa + bs2*ai[row]; 654 rmax = imax[row]; 655 nrow = ailen[row]; 656 low = 0; 657 high = nrow; 658 for (l=0; l<n; l++) { /* loop over added columns */ 659 if (in[l] < 0) continue; 660 col = in[l]; 661 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); 662 if (col < row) { 663 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */ 664 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)"); 665 } 666 if (roworiented) value = v + k*(stepval+bs)*bs + l*bs; 667 else value = v + l*(stepval+bs)*bs + k*bs; 668 669 if (col <= lastcol) low = 0; 670 else high = nrow; 671 672 lastcol = col; 673 while (high-low > 7) { 674 t = (low+high)/2; 675 if (rp[t] > col) high = t; 676 else low = t; 677 } 678 for (i=low; i<high; i++) { 679 if (rp[i] > col) break; 680 if (rp[i] == col) { 681 bap = ap + bs2*i; 682 if (roworiented) { 683 if (is == ADD_VALUES) { 684 for (ii=0; ii<bs; ii++,value+=stepval) { 685 for (jj=ii; jj<bs2; jj+=bs) { 686 bap[jj] += *value++; 687 } 688 } 689 } else { 690 for (ii=0; ii<bs; ii++,value+=stepval) { 691 for (jj=ii; jj<bs2; jj+=bs) { 692 bap[jj] = *value++; 693 } 694 } 695 } 696 } else { 697 if (is == ADD_VALUES) { 698 for (ii=0; ii<bs; ii++,value+=stepval) { 699 for (jj=0; jj<bs; jj++) { 700 *bap++ += *value++; 701 } 702 } 703 } else { 704 for (ii=0; ii<bs; ii++,value+=stepval) { 705 for (jj=0; jj<bs; jj++) { 706 *bap++ = *value++; 707 } 708 } 709 } 710 } 711 goto noinsert2; 712 } 713 } 714 if (nonew == 1) goto noinsert2; 715 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); 716 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 717 N = nrow++ - 1; high++; 718 /* shift up all the later entries in this row */ 719 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 720 PetscCall(PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1))); 721 PetscCall(PetscArrayzero(ap+bs2*i,bs2)); 722 rp[i] = col; 723 bap = ap + bs2*i; 724 if (roworiented) { 725 for (ii=0; ii<bs; ii++,value+=stepval) { 726 for (jj=ii; jj<bs2; jj+=bs) { 727 bap[jj] = *value++; 728 } 729 } 730 } else { 731 for (ii=0; ii<bs; ii++,value+=stepval) { 732 for (jj=0; jj<bs; jj++) { 733 *bap++ = *value++; 734 } 735 } 736 } 737 noinsert2:; 738 low = i; 739 } 740 ailen[row] = nrow; 741 } 742 PetscFunctionReturn(0); 743 } 744 745 /* 746 This is not yet used 747 */ 748 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 749 { 750 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 751 const PetscInt *ai = a->i, *aj = a->j,*cols; 752 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 753 PetscBool flag; 754 755 PetscFunctionBegin; 756 PetscCall(PetscMalloc1(m,&ns)); 757 while (i < m) { 758 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 759 /* Limits the number of elements in a node to 'a->inode.limit' */ 760 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 761 nzy = ai[j+1] - ai[j]; 762 if (nzy != (nzx - j + i)) break; 763 PetscCall(PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag)); 764 if (!flag) break; 765 } 766 ns[node_count++] = blk_size; 767 768 i = j; 769 } 770 if (!a->inode.size && m && node_count > .9*m) { 771 PetscCall(PetscFree(ns)); 772 PetscCall(PetscInfo(A,"Found %" PetscInt_FMT " nodes out of %" PetscInt_FMT " rows. Not using Inode routines\n",node_count,m)); 773 } else { 774 a->inode.node_count = node_count; 775 776 PetscCall(PetscMalloc1(node_count,&a->inode.size)); 777 PetscCall(PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt))); 778 PetscCall(PetscArraycpy(a->inode.size,ns,node_count)); 779 PetscCall(PetscFree(ns)); 780 PetscCall(PetscInfo(A,"Found %" PetscInt_FMT " nodes of %" PetscInt_FMT ". Limit used: %" PetscInt_FMT ". Using Inode routines\n",node_count,m,a->inode.limit)); 781 782 /* count collections of adjacent columns in each inode */ 783 row = 0; 784 cnt = 0; 785 for (i=0; i<node_count; i++) { 786 cols = aj + ai[row] + a->inode.size[i]; 787 nz = ai[row+1] - ai[row] - a->inode.size[i]; 788 for (j=1; j<nz; j++) { 789 if (cols[j] != cols[j-1]+1) cnt++; 790 } 791 cnt++; 792 row += a->inode.size[i]; 793 } 794 PetscCall(PetscMalloc1(2*cnt,&counts)); 795 cnt = 0; 796 row = 0; 797 for (i=0; i<node_count; i++) { 798 cols = aj + ai[row] + a->inode.size[i]; 799 counts[2*cnt] = cols[0]; 800 nz = ai[row+1] - ai[row] - a->inode.size[i]; 801 cnt2 = 1; 802 for (j=1; j<nz; j++) { 803 if (cols[j] != cols[j-1]+1) { 804 counts[2*(cnt++)+1] = cnt2; 805 counts[2*cnt] = cols[j]; 806 cnt2 = 1; 807 } else cnt2++; 808 } 809 counts[2*(cnt++)+1] = cnt2; 810 row += a->inode.size[i]; 811 } 812 PetscCall(PetscIntView(2*cnt,counts,NULL)); 813 } 814 PetscFunctionReturn(0); 815 } 816 817 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 818 { 819 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 820 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 821 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 822 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 823 MatScalar *aa = a->a,*ap; 824 825 PetscFunctionBegin; 826 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 827 828 if (m) rmax = ailen[0]; 829 for (i=1; i<mbs; i++) { 830 /* move each row back by the amount of empty slots (fshift) before it*/ 831 fshift += imax[i-1] - ailen[i-1]; 832 rmax = PetscMax(rmax,ailen[i]); 833 if (fshift) { 834 ip = aj + ai[i]; 835 ap = aa + bs2*ai[i]; 836 N = ailen[i]; 837 PetscCall(PetscArraymove(ip-fshift,ip,N)); 838 PetscCall(PetscArraymove(ap-bs2*fshift,ap,bs2*N)); 839 } 840 ai[i] = ai[i-1] + ailen[i-1]; 841 } 842 if (mbs) { 843 fshift += imax[mbs-1] - ailen[mbs-1]; 844 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 845 } 846 /* reset ilen and imax for each row */ 847 for (i=0; i<mbs; i++) { 848 ailen[i] = imax[i] = ai[i+1] - ai[i]; 849 } 850 a->nz = ai[mbs]; 851 852 /* diagonals may have moved, reset it */ 853 if (a->diag) { 854 PetscCall(PetscArraycpy(a->diag,ai,mbs)); 855 } 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 }; 1516 1517 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1518 { 1519 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1520 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1521 1522 PetscFunctionBegin; 1523 PetscCheck(aij->nonew == 1,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1524 1525 /* allocate space for values if not already there */ 1526 if (!aij->saved_values) { 1527 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 1528 } 1529 1530 /* copy values over */ 1531 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 1532 PetscFunctionReturn(0); 1533 } 1534 1535 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1536 { 1537 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1538 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1539 1540 PetscFunctionBegin; 1541 PetscCheck(aij->nonew == 1,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1542 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1543 1544 /* copy values over */ 1545 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 1546 PetscFunctionReturn(0); 1547 } 1548 1549 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1550 { 1551 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1552 PetscInt i,mbs,nbs,bs2; 1553 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1554 1555 PetscFunctionBegin; 1556 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1557 1558 PetscCall(MatSetBlockSize(B,PetscAbs(bs))); 1559 PetscCall(PetscLayoutSetUp(B->rmap)); 1560 PetscCall(PetscLayoutSetUp(B->cmap)); 1561 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); 1562 PetscCall(PetscLayoutGetBlockSize(B->rmap,&bs)); 1563 1564 B->preallocated = PETSC_TRUE; 1565 1566 mbs = B->rmap->N/bs; 1567 nbs = B->cmap->n/bs; 1568 bs2 = bs*bs; 1569 1570 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"); 1571 1572 if (nz == MAT_SKIP_ALLOCATION) { 1573 skipallocation = PETSC_TRUE; 1574 nz = 0; 1575 } 1576 1577 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1578 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 1579 if (nnz) { 1580 for (i=0; i<mbs; i++) { 1581 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]); 1582 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); 1583 } 1584 } 1585 1586 B->ops->mult = MatMult_SeqSBAIJ_N; 1587 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1588 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1589 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1590 1591 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL)); 1592 if (!flg) { 1593 switch (bs) { 1594 case 1: 1595 B->ops->mult = MatMult_SeqSBAIJ_1; 1596 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1597 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1598 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1599 break; 1600 case 2: 1601 B->ops->mult = MatMult_SeqSBAIJ_2; 1602 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1603 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1604 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1605 break; 1606 case 3: 1607 B->ops->mult = MatMult_SeqSBAIJ_3; 1608 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1609 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1610 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1611 break; 1612 case 4: 1613 B->ops->mult = MatMult_SeqSBAIJ_4; 1614 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1615 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1616 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1617 break; 1618 case 5: 1619 B->ops->mult = MatMult_SeqSBAIJ_5; 1620 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1621 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1622 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1623 break; 1624 case 6: 1625 B->ops->mult = MatMult_SeqSBAIJ_6; 1626 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1627 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1628 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1629 break; 1630 case 7: 1631 B->ops->mult = MatMult_SeqSBAIJ_7; 1632 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1633 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1634 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1635 break; 1636 } 1637 } 1638 1639 b->mbs = mbs; 1640 b->nbs = nbs; 1641 if (!skipallocation) { 1642 if (!b->imax) { 1643 PetscCall(PetscMalloc2(mbs,&b->imax,mbs,&b->ilen)); 1644 1645 b->free_imax_ilen = PETSC_TRUE; 1646 1647 PetscCall(PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt))); 1648 } 1649 if (!nnz) { 1650 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1651 else if (nz <= 0) nz = 1; 1652 nz = PetscMin(nbs,nz); 1653 for (i=0; i<mbs; i++) b->imax[i] = nz; 1654 PetscCall(PetscIntMultError(nz,mbs,&nz)); 1655 } else { 1656 PetscInt64 nz64 = 0; 1657 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 1658 PetscCall(PetscIntCast(nz64,&nz)); 1659 } 1660 /* b->ilen will count nonzeros in each block row so far. */ 1661 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1662 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1663 1664 /* allocate the matrix space */ 1665 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 1666 PetscCall(PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i)); 1667 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)))); 1668 PetscCall(PetscArrayzero(b->a,nz*bs2)); 1669 PetscCall(PetscArrayzero(b->j,nz)); 1670 1671 b->singlemalloc = PETSC_TRUE; 1672 1673 /* pointer to beginning of each row */ 1674 b->i[0] = 0; 1675 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1676 1677 b->free_a = PETSC_TRUE; 1678 b->free_ij = PETSC_TRUE; 1679 } else { 1680 b->free_a = PETSC_FALSE; 1681 b->free_ij = PETSC_FALSE; 1682 } 1683 1684 b->bs2 = bs2; 1685 b->nz = 0; 1686 b->maxnz = nz; 1687 b->inew = NULL; 1688 b->jnew = NULL; 1689 b->anew = NULL; 1690 b->a2anew = NULL; 1691 b->permute = PETSC_FALSE; 1692 1693 B->was_assembled = PETSC_FALSE; 1694 B->assembled = PETSC_FALSE; 1695 if (realalloc) PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 1696 PetscFunctionReturn(0); 1697 } 1698 1699 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1700 { 1701 PetscInt i,j,m,nz,anz, nz_max=0,*nnz; 1702 PetscScalar *values=NULL; 1703 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1704 1705 PetscFunctionBegin; 1706 PetscCheck(bs >= 1,PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %" PetscInt_FMT,bs); 1707 PetscCall(PetscLayoutSetBlockSize(B->rmap,bs)); 1708 PetscCall(PetscLayoutSetBlockSize(B->cmap,bs)); 1709 PetscCall(PetscLayoutSetUp(B->rmap)); 1710 PetscCall(PetscLayoutSetUp(B->cmap)); 1711 PetscCall(PetscLayoutGetBlockSize(B->rmap,&bs)); 1712 m = B->rmap->n/bs; 1713 1714 PetscCheck(!ii[0],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %" PetscInt_FMT,ii[0]); 1715 PetscCall(PetscMalloc1(m+1,&nnz)); 1716 for (i=0; i<m; i++) { 1717 nz = ii[i+1] - ii[i]; 1718 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT,i,nz); 1719 anz = 0; 1720 for (j=0; j<nz; j++) { 1721 /* count only values on the diagonal or above */ 1722 if (jj[ii[i] + j] >= i) { 1723 anz = nz - j; 1724 break; 1725 } 1726 } 1727 nz_max = PetscMax(nz_max,anz); 1728 nnz[i] = anz; 1729 } 1730 PetscCall(MatSeqSBAIJSetPreallocation(B,bs,0,nnz)); 1731 PetscCall(PetscFree(nnz)); 1732 1733 values = (PetscScalar*)V; 1734 if (!values) { 1735 PetscCall(PetscCalloc1(bs*bs*nz_max,&values)); 1736 } 1737 for (i=0; i<m; i++) { 1738 PetscInt ncols = ii[i+1] - ii[i]; 1739 const PetscInt *icols = jj + ii[i]; 1740 if (!roworiented || bs == 1) { 1741 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1742 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES)); 1743 } else { 1744 for (j=0; j<ncols; j++) { 1745 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1746 PetscCall(MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES)); 1747 } 1748 } 1749 } 1750 if (!V) PetscCall(PetscFree(values)); 1751 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1752 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1753 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 1754 PetscFunctionReturn(0); 1755 } 1756 1757 /* 1758 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1759 */ 1760 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1761 { 1762 PetscBool flg = PETSC_FALSE; 1763 PetscInt bs = B->rmap->bs; 1764 1765 PetscFunctionBegin; 1766 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL)); 1767 if (flg) bs = 8; 1768 1769 if (!natural) { 1770 switch (bs) { 1771 case 1: 1772 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1773 break; 1774 case 2: 1775 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1776 break; 1777 case 3: 1778 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1779 break; 1780 case 4: 1781 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1782 break; 1783 case 5: 1784 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1785 break; 1786 case 6: 1787 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1788 break; 1789 case 7: 1790 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1791 break; 1792 default: 1793 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1794 break; 1795 } 1796 } else { 1797 switch (bs) { 1798 case 1: 1799 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1800 break; 1801 case 2: 1802 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1803 break; 1804 case 3: 1805 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1806 break; 1807 case 4: 1808 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1809 break; 1810 case 5: 1811 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1812 break; 1813 case 6: 1814 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1815 break; 1816 case 7: 1817 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1818 break; 1819 default: 1820 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1821 break; 1822 } 1823 } 1824 PetscFunctionReturn(0); 1825 } 1826 1827 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*); 1828 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*); 1829 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type) 1830 { 1831 PetscFunctionBegin; 1832 *type = MATSOLVERPETSC; 1833 PetscFunctionReturn(0); 1834 } 1835 1836 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1837 { 1838 PetscInt n = A->rmap->n; 1839 1840 PetscFunctionBegin; 1841 #if defined(PETSC_USE_COMPLEX) 1842 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"); 1843 #endif 1844 1845 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 1846 PetscCall(MatSetSizes(*B,n,n,n,n)); 1847 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1848 PetscCall(MatSetType(*B,MATSEQSBAIJ)); 1849 PetscCall(MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL)); 1850 1851 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1852 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1853 PetscCall(PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY])); 1854 PetscCall(PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC])); 1855 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1856 1857 (*B)->factortype = ftype; 1858 (*B)->canuseordering = PETSC_TRUE; 1859 PetscCall(PetscFree((*B)->solvertype)); 1860 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype)); 1861 PetscCall(PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc)); 1862 PetscFunctionReturn(0); 1863 } 1864 1865 /*@C 1866 MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored 1867 1868 Not Collective 1869 1870 Input Parameter: 1871 . mat - a MATSEQSBAIJ matrix 1872 1873 Output Parameter: 1874 . array - pointer to the data 1875 1876 Level: intermediate 1877 1878 .seealso: `MatSeqSBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1879 @*/ 1880 PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array) 1881 { 1882 PetscFunctionBegin; 1883 PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array)); 1884 PetscFunctionReturn(0); 1885 } 1886 1887 /*@C 1888 MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray() 1889 1890 Not Collective 1891 1892 Input Parameters: 1893 + mat - a MATSEQSBAIJ matrix 1894 - array - pointer to the data 1895 1896 Level: intermediate 1897 1898 .seealso: `MatSeqSBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()` 1899 @*/ 1900 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array) 1901 { 1902 PetscFunctionBegin; 1903 PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array)); 1904 PetscFunctionReturn(0); 1905 } 1906 1907 /*MC 1908 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1909 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1910 1911 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1912 can call MatSetOption(Mat, MAT_HERMITIAN). 1913 1914 Options Database Keys: 1915 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1916 1917 Notes: 1918 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1919 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1920 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. 1921 1922 The number of rows in the matrix must be less than or equal to the number of columns 1923 1924 Level: beginner 1925 1926 .seealso: `MatCreateSeqSBAIJ()`, `MatType`, `MATMPISBAIJ` 1927 M*/ 1928 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1929 { 1930 Mat_SeqSBAIJ *b; 1931 PetscMPIInt size; 1932 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1933 1934 PetscFunctionBegin; 1935 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 1936 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1937 1938 PetscCall(PetscNewLog(B,&b)); 1939 B->data = (void*)b; 1940 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 1941 1942 B->ops->destroy = MatDestroy_SeqSBAIJ; 1943 B->ops->view = MatView_SeqSBAIJ; 1944 b->row = NULL; 1945 b->icol = NULL; 1946 b->reallocs = 0; 1947 b->saved_values = NULL; 1948 b->inode.limit = 5; 1949 b->inode.max_limit = 5; 1950 1951 b->roworiented = PETSC_TRUE; 1952 b->nonew = 0; 1953 b->diag = NULL; 1954 b->solve_work = NULL; 1955 b->mult_work = NULL; 1956 B->spptr = NULL; 1957 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1958 b->keepnonzeropattern = PETSC_FALSE; 1959 1960 b->inew = NULL; 1961 b->jnew = NULL; 1962 b->anew = NULL; 1963 b->a2anew = NULL; 1964 b->permute = PETSC_FALSE; 1965 1966 b->ignore_ltriangular = PETSC_TRUE; 1967 1968 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL)); 1969 1970 b->getrow_utriangular = PETSC_FALSE; 1971 1972 PetscCall(PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL)); 1973 1974 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ)); 1975 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ)); 1976 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ)); 1977 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ)); 1978 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ)); 1979 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ)); 1980 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ)); 1981 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ)); 1982 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ)); 1983 #if defined(PETSC_HAVE_ELEMENTAL) 1984 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental)); 1985 #endif 1986 #if defined(PETSC_HAVE_SCALAPACK) 1987 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK)); 1988 #endif 1989 1990 B->symmetric = PETSC_TRUE; 1991 B->structurally_symmetric = PETSC_TRUE; 1992 B->symmetric_set = PETSC_TRUE; 1993 B->structurally_symmetric_set = PETSC_TRUE; 1994 B->symmetric_eternal = PETSC_TRUE; 1995 #if defined(PETSC_USE_COMPLEX) 1996 B->hermitian = PETSC_FALSE; 1997 B->hermitian_set = PETSC_FALSE; 1998 #else 1999 B->hermitian = PETSC_TRUE; 2000 B->hermitian_set = PETSC_TRUE; 2001 #endif 2002 2003 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ)); 2004 2005 PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat"); 2006 PetscCall(PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL)); 2007 if (no_unroll) { 2008 PetscCall(PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n")); 2009 } 2010 PetscCall(PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL)); 2011 if (no_inode) PetscCall(PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n")); 2012 PetscCall(PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL)); 2013 PetscOptionsEnd(); 2014 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 2015 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2016 PetscFunctionReturn(0); 2017 } 2018 2019 /*@C 2020 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2021 compressed row) format. For good matrix assembly performance the 2022 user should preallocate the matrix storage by setting the parameter nz 2023 (or the array nnz). By setting these parameters accurately, performance 2024 during matrix assembly can be increased by more than a factor of 50. 2025 2026 Collective on Mat 2027 2028 Input Parameters: 2029 + B - the symmetric matrix 2030 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2031 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2032 . nz - number of block nonzeros per block row (same for all rows) 2033 - nnz - array containing the number of block nonzeros in the upper triangular plus 2034 diagonal portion of each block (possibly different for each block row) or NULL 2035 2036 Options Database Keys: 2037 + -mat_no_unroll - uses code that does not unroll the loops in the 2038 block calculations (much slower) 2039 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2040 2041 Level: intermediate 2042 2043 Notes: 2044 Specify the preallocated storage with either nz or nnz (not both). 2045 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2046 allocation. See Users-Manual: ch_mat for details. 2047 2048 You can call MatGetInfo() to get information on how effective the preallocation was; 2049 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2050 You can also run with the option -info and look for messages with the string 2051 malloc in them to see if additional memory allocation was needed. 2052 2053 If the nnz parameter is given then the nz parameter is ignored 2054 2055 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2056 @*/ 2057 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2058 { 2059 PetscFunctionBegin; 2060 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2061 PetscValidType(B,1); 2062 PetscValidLogicalCollectiveInt(B,bs,2); 2063 PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz)); 2064 PetscFunctionReturn(0); 2065 } 2066 2067 /*@C 2068 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values 2069 2070 Input Parameters: 2071 + B - the matrix 2072 . bs - size of block, the blocks are ALWAYS square. 2073 . i - the indices into j for the start of each local row (starts with zero) 2074 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2075 - v - optional values in the matrix 2076 2077 Level: advanced 2078 2079 Notes: 2080 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2081 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2082 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2083 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2084 block column and the second index is over columns within a block. 2085 2086 Any entries below the diagonal are ignored 2087 2088 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 2089 and usually the numerical values as well 2090 2091 .seealso: `MatCreate()`, `MatCreateSeqSBAIJ()`, `MatSetValuesBlocked()`, `MatSeqSBAIJSetPreallocation()`, `MATSEQSBAIJ` 2092 @*/ 2093 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2094 { 2095 PetscFunctionBegin; 2096 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2097 PetscValidType(B,1); 2098 PetscValidLogicalCollectiveInt(B,bs,2); 2099 PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v)); 2100 PetscFunctionReturn(0); 2101 } 2102 2103 /*@C 2104 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2105 compressed row) format. For good matrix assembly performance the 2106 user should preallocate the matrix storage by setting the parameter nz 2107 (or the array nnz). By setting these parameters accurately, performance 2108 during matrix assembly can be increased by more than a factor of 50. 2109 2110 Collective 2111 2112 Input Parameters: 2113 + comm - MPI communicator, set to PETSC_COMM_SELF 2114 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2115 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2116 . m - number of rows, or number of columns 2117 . nz - number of block nonzeros per block row (same for all rows) 2118 - nnz - array containing the number of block nonzeros in the upper triangular plus 2119 diagonal portion of each block (possibly different for each block row) or NULL 2120 2121 Output Parameter: 2122 . A - the symmetric matrix 2123 2124 Options Database Keys: 2125 + -mat_no_unroll - uses code that does not unroll the loops in the 2126 block calculations (much slower) 2127 - -mat_block_size - size of the blocks to use 2128 2129 Level: intermediate 2130 2131 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2132 MatXXXXSetPreallocation() paradigm instead of this routine directly. 2133 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2134 2135 Notes: 2136 The number of rows and columns must be divisible by blocksize. 2137 This matrix type does not support complex Hermitian operation. 2138 2139 Specify the preallocated storage with either nz or nnz (not both). 2140 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2141 allocation. See Users-Manual: ch_mat for details. 2142 2143 If the nnz parameter is given then the nz parameter is ignored 2144 2145 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateSBAIJ()` 2146 @*/ 2147 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2148 { 2149 PetscFunctionBegin; 2150 PetscCall(MatCreate(comm,A)); 2151 PetscCall(MatSetSizes(*A,m,n,m,n)); 2152 PetscCall(MatSetType(*A,MATSEQSBAIJ)); 2153 PetscCall(MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz)); 2154 PetscFunctionReturn(0); 2155 } 2156 2157 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2158 { 2159 Mat C; 2160 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2161 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2162 2163 PetscFunctionBegin; 2164 PetscCheck(a->i[mbs] == nz,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2165 2166 *B = NULL; 2167 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2168 PetscCall(MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n)); 2169 PetscCall(MatSetBlockSizesFromMats(C,A,A)); 2170 PetscCall(MatSetType(C,MATSEQSBAIJ)); 2171 c = (Mat_SeqSBAIJ*)C->data; 2172 2173 C->preallocated = PETSC_TRUE; 2174 C->factortype = A->factortype; 2175 c->row = NULL; 2176 c->icol = NULL; 2177 c->saved_values = NULL; 2178 c->keepnonzeropattern = a->keepnonzeropattern; 2179 C->assembled = PETSC_TRUE; 2180 2181 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 2182 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 2183 c->bs2 = a->bs2; 2184 c->mbs = a->mbs; 2185 c->nbs = a->nbs; 2186 2187 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2188 c->imax = a->imax; 2189 c->ilen = a->ilen; 2190 c->free_imax_ilen = PETSC_FALSE; 2191 } else { 2192 PetscCall(PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen)); 2193 PetscCall(PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt))); 2194 for (i=0; i<mbs; i++) { 2195 c->imax[i] = a->imax[i]; 2196 c->ilen[i] = a->ilen[i]; 2197 } 2198 c->free_imax_ilen = PETSC_TRUE; 2199 } 2200 2201 /* allocate the matrix space */ 2202 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2203 PetscCall(PetscMalloc1(bs2*nz,&c->a)); 2204 PetscCall(PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar))); 2205 c->i = a->i; 2206 c->j = a->j; 2207 c->singlemalloc = PETSC_FALSE; 2208 c->free_a = PETSC_TRUE; 2209 c->free_ij = PETSC_FALSE; 2210 c->parent = A; 2211 PetscCall(PetscObjectReference((PetscObject)A)); 2212 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 2213 PetscCall(MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 2214 } else { 2215 PetscCall(PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i)); 2216 PetscCall(PetscArraycpy(c->i,a->i,mbs+1)); 2217 PetscCall(PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)))); 2218 c->singlemalloc = PETSC_TRUE; 2219 c->free_a = PETSC_TRUE; 2220 c->free_ij = PETSC_TRUE; 2221 } 2222 if (mbs > 0) { 2223 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2224 PetscCall(PetscArraycpy(c->j,a->j,nz)); 2225 } 2226 if (cpvalues == MAT_COPY_VALUES) { 2227 PetscCall(PetscArraycpy(c->a,a->a,bs2*nz)); 2228 } else { 2229 PetscCall(PetscArrayzero(c->a,bs2*nz)); 2230 } 2231 if (a->jshort) { 2232 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2233 /* if the parent matrix is reassembled, this child matrix will never notice */ 2234 PetscCall(PetscMalloc1(nz,&c->jshort)); 2235 PetscCall(PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short))); 2236 PetscCall(PetscArraycpy(c->jshort,a->jshort,nz)); 2237 2238 c->free_jshort = PETSC_TRUE; 2239 } 2240 } 2241 2242 c->roworiented = a->roworiented; 2243 c->nonew = a->nonew; 2244 2245 if (a->diag) { 2246 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2247 c->diag = a->diag; 2248 c->free_diag = PETSC_FALSE; 2249 } else { 2250 PetscCall(PetscMalloc1(mbs,&c->diag)); 2251 PetscCall(PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt))); 2252 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2253 c->free_diag = PETSC_TRUE; 2254 } 2255 } 2256 c->nz = a->nz; 2257 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2258 c->solve_work = NULL; 2259 c->mult_work = NULL; 2260 2261 *B = C; 2262 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 2263 PetscFunctionReturn(0); 2264 } 2265 2266 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2267 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2268 2269 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer) 2270 { 2271 PetscBool isbinary; 2272 2273 PetscFunctionBegin; 2274 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 2275 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); 2276 PetscCall(MatLoad_SeqSBAIJ_Binary(mat,viewer)); 2277 PetscFunctionReturn(0); 2278 } 2279 2280 /*@ 2281 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2282 (upper triangular entries in CSR format) provided by the user. 2283 2284 Collective 2285 2286 Input Parameters: 2287 + comm - must be an MPI communicator of size 1 2288 . bs - size of block 2289 . m - number of rows 2290 . n - number of columns 2291 . 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 2292 . j - column indices 2293 - a - matrix values 2294 2295 Output Parameter: 2296 . mat - the matrix 2297 2298 Level: advanced 2299 2300 Notes: 2301 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2302 once the matrix is destroyed 2303 2304 You cannot set new nonzero locations into this matrix, that will generate an error. 2305 2306 The i and j indices are 0 based 2307 2308 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 2309 it is the regular CSR format excluding the lower triangular elements. 2310 2311 .seealso: `MatCreate()`, `MatCreateSBAIJ()`, `MatCreateSeqSBAIJ()` 2312 2313 @*/ 2314 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 2315 { 2316 PetscInt ii; 2317 Mat_SeqSBAIJ *sbaij; 2318 2319 PetscFunctionBegin; 2320 PetscCheck(bs == 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %" PetscInt_FMT " > 1 is not supported yet",bs); 2321 PetscCheck(m == 0 || i[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2322 2323 PetscCall(MatCreate(comm,mat)); 2324 PetscCall(MatSetSizes(*mat,m,n,m,n)); 2325 PetscCall(MatSetType(*mat,MATSEQSBAIJ)); 2326 PetscCall(MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL)); 2327 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2328 PetscCall(PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen)); 2329 PetscCall(PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt))); 2330 2331 sbaij->i = i; 2332 sbaij->j = j; 2333 sbaij->a = a; 2334 2335 sbaij->singlemalloc = PETSC_FALSE; 2336 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2337 sbaij->free_a = PETSC_FALSE; 2338 sbaij->free_ij = PETSC_FALSE; 2339 sbaij->free_imax_ilen = PETSC_TRUE; 2340 2341 for (ii=0; ii<m; ii++) { 2342 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2343 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]); 2344 } 2345 if (PetscDefined(USE_DEBUG)) { 2346 for (ii=0; ii<sbaij->i[m]; ii++) { 2347 PetscCheck(j[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 2348 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]); 2349 } 2350 } 2351 2352 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 2353 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 2354 PetscFunctionReturn(0); 2355 } 2356 2357 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2358 { 2359 PetscMPIInt size; 2360 2361 PetscFunctionBegin; 2362 PetscCallMPI(MPI_Comm_size(comm,&size)); 2363 if (size == 1 && scall == MAT_REUSE_MATRIX) { 2364 PetscCall(MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN)); 2365 } else { 2366 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat)); 2367 } 2368 PetscFunctionReturn(0); 2369 } 2370