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