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