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