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