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