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