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