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