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