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