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