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