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