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