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