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_SeqSBAIJ(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 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 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1739 PetscFunctionReturn(0); 1740 } 1741 1742 #undef __FUNCT__ 1743 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ" 1744 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1745 { 1746 PetscInt i,j,m,nz,nz_max=0,*nnz; 1747 PetscScalar *values=0; 1748 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1749 PetscErrorCode ierr; 1750 PetscFunctionBegin; 1751 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs); 1752 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1753 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1754 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1755 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1756 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1757 m = B->rmap->n/bs; 1758 1759 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 1760 ierr = PetscMalloc1(m+1,&nnz);CHKERRQ(ierr); 1761 for (i=0; i<m; i++) { 1762 nz = ii[i+1] - ii[i]; 1763 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz); 1764 nz_max = PetscMax(nz_max,nz); 1765 nnz[i] = nz; 1766 } 1767 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1768 ierr = PetscFree(nnz);CHKERRQ(ierr); 1769 1770 values = (PetscScalar*)V; 1771 if (!values) { 1772 ierr = PetscCalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr); 1773 } 1774 for (i=0; i<m; i++) { 1775 PetscInt ncols = ii[i+1] - ii[i]; 1776 const PetscInt *icols = jj + ii[i]; 1777 if (!roworiented || bs == 1) { 1778 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1779 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 1780 } else { 1781 for (j=0; j<ncols; j++) { 1782 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1783 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 1784 } 1785 } 1786 } 1787 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 1788 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1789 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1790 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1791 PetscFunctionReturn(0); 1792 } 1793 1794 /* 1795 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1796 */ 1797 #undef __FUNCT__ 1798 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1799 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1800 { 1801 PetscErrorCode ierr; 1802 PetscBool flg = PETSC_FALSE; 1803 PetscInt bs = B->rmap->bs; 1804 1805 PetscFunctionBegin; 1806 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1807 if (flg) bs = 8; 1808 1809 if (!natural) { 1810 switch (bs) { 1811 case 1: 1812 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1813 break; 1814 case 2: 1815 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1816 break; 1817 case 3: 1818 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1819 break; 1820 case 4: 1821 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1822 break; 1823 case 5: 1824 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1825 break; 1826 case 6: 1827 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1828 break; 1829 case 7: 1830 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1831 break; 1832 default: 1833 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1834 break; 1835 } 1836 } else { 1837 switch (bs) { 1838 case 1: 1839 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1840 break; 1841 case 2: 1842 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1843 break; 1844 case 3: 1845 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1846 break; 1847 case 4: 1848 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1849 break; 1850 case 5: 1851 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1852 break; 1853 case 6: 1854 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1855 break; 1856 case 7: 1857 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1858 break; 1859 default: 1860 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1861 break; 1862 } 1863 } 1864 PetscFunctionReturn(0); 1865 } 1866 1867 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1868 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1869 1870 #undef __FUNCT__ 1871 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1872 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1873 { 1874 PetscInt n = A->rmap->n; 1875 PetscErrorCode ierr; 1876 1877 PetscFunctionBegin; 1878 #if defined(PETSC_USE_COMPLEX) 1879 if (A->hermitian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 1880 #endif 1881 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1882 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1883 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1884 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1885 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1886 1887 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1888 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1889 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1890 1891 (*B)->factortype = ftype; 1892 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 1893 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 1894 PetscFunctionReturn(0); 1895 } 1896 1897 /*MC 1898 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1899 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1900 1901 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1902 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1903 1904 Options Database Keys: 1905 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1906 1907 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1908 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1909 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. 1910 1911 1912 Level: beginner 1913 1914 .seealso: MatCreateSeqSBAIJ 1915 M*/ 1916 1917 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1918 1919 #undef __FUNCT__ 1920 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1921 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1922 { 1923 Mat_SeqSBAIJ *b; 1924 PetscErrorCode ierr; 1925 PetscMPIInt size; 1926 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1927 1928 PetscFunctionBegin; 1929 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 1930 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1931 1932 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1933 B->data = (void*)b; 1934 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1935 1936 B->ops->destroy = MatDestroy_SeqSBAIJ; 1937 B->ops->view = MatView_SeqSBAIJ; 1938 b->row = 0; 1939 b->icol = 0; 1940 b->reallocs = 0; 1941 b->saved_values = 0; 1942 b->inode.limit = 5; 1943 b->inode.max_limit = 5; 1944 1945 b->roworiented = PETSC_TRUE; 1946 b->nonew = 0; 1947 b->diag = 0; 1948 b->solve_work = 0; 1949 b->mult_work = 0; 1950 B->spptr = 0; 1951 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1952 b->keepnonzeropattern = PETSC_FALSE; 1953 1954 b->inew = 0; 1955 b->jnew = 0; 1956 b->anew = 0; 1957 b->a2anew = 0; 1958 b->permute = PETSC_FALSE; 1959 1960 b->ignore_ltriangular = PETSC_TRUE; 1961 1962 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1963 1964 b->getrow_utriangular = PETSC_FALSE; 1965 1966 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 1967 1968 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1969 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1970 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1971 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1972 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1973 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1974 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 1975 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C",MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1976 #if defined(PETSC_HAVE_ELEMENTAL) 1977 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr); 1978 #endif 1979 1980 B->symmetric = PETSC_TRUE; 1981 B->structurally_symmetric = PETSC_TRUE; 1982 B->symmetric_set = PETSC_TRUE; 1983 B->structurally_symmetric_set = PETSC_TRUE; 1984 1985 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1986 1987 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1988 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 1989 if (no_unroll) { 1990 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 1991 } 1992 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 1993 if (no_inode) { 1994 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 1995 } 1996 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 1997 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1998 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1999 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2000 PetscFunctionReturn(0); 2001 } 2002 2003 #undef __FUNCT__ 2004 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 2005 /*@C 2006 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2007 compressed row) format. For good matrix assembly performance the 2008 user should preallocate the matrix storage by setting the parameter nz 2009 (or the array nnz). By setting these parameters accurately, performance 2010 during matrix assembly can be increased by more than a factor of 50. 2011 2012 Collective on Mat 2013 2014 Input Parameters: 2015 + B - the symmetric matrix 2016 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2017 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2018 . nz - number of block nonzeros per block row (same for all rows) 2019 - nnz - array containing the number of block nonzeros in the upper triangular plus 2020 diagonal portion of each block (possibly different for each block row) or NULL 2021 2022 Options Database Keys: 2023 . -mat_no_unroll - uses code that does not unroll the loops in the 2024 block calculations (much slower) 2025 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2026 2027 Level: intermediate 2028 2029 Notes: 2030 Specify the preallocated storage with either nz or nnz (not both). 2031 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2032 allocation. See Users-Manual: ch_mat for details. 2033 2034 You can call MatGetInfo() to get information on how effective the preallocation was; 2035 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2036 You can also run with the option -info and look for messages with the string 2037 malloc in them to see if additional memory allocation was needed. 2038 2039 If the nnz parameter is given then the nz parameter is ignored 2040 2041 2042 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2043 @*/ 2044 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2045 { 2046 PetscErrorCode ierr; 2047 2048 PetscFunctionBegin; 2049 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2050 PetscValidType(B,1); 2051 PetscValidLogicalCollectiveInt(B,bs,2); 2052 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2053 PetscFunctionReturn(0); 2054 } 2055 2056 #undef __FUNCT__ 2057 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR" 2058 /*@C 2059 MatSeqSBAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in symmetric block AIJ format. 2060 2061 Input Parameters: 2062 + B - the matrix 2063 . bs - size of block, the blocks are ALWAYS square. 2064 . i - the indices into j for the start of each local row (starts with zero) 2065 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2066 - v - optional values in the matrix 2067 2068 Level: developer 2069 2070 Notes: 2071 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2072 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2073 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2074 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2075 block column and the second index is over columns within a block. 2076 2077 .keywords: matrix, block, aij, compressed row, sparse 2078 2079 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2080 @*/ 2081 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2082 { 2083 PetscErrorCode ierr; 2084 2085 PetscFunctionBegin; 2086 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2087 PetscValidType(B,1); 2088 PetscValidLogicalCollectiveInt(B,bs,2); 2089 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2090 PetscFunctionReturn(0); 2091 } 2092 2093 #undef __FUNCT__ 2094 #define __FUNCT__ "MatCreateSeqSBAIJ" 2095 /*@C 2096 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2097 compressed row) format. For good matrix assembly performance the 2098 user should preallocate the matrix storage by setting the parameter nz 2099 (or the array nnz). By setting these parameters accurately, performance 2100 during matrix assembly can be increased by more than a factor of 50. 2101 2102 Collective on MPI_Comm 2103 2104 Input Parameters: 2105 + comm - MPI communicator, set to PETSC_COMM_SELF 2106 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2107 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2108 . m - number of rows, or number of columns 2109 . nz - number of block nonzeros per block row (same for all rows) 2110 - nnz - array containing the number of block nonzeros in the upper triangular plus 2111 diagonal portion of each block (possibly different for each block row) or NULL 2112 2113 Output Parameter: 2114 . A - the symmetric matrix 2115 2116 Options Database Keys: 2117 . -mat_no_unroll - uses code that does not unroll the loops in the 2118 block calculations (much slower) 2119 . -mat_block_size - size of the blocks to use 2120 2121 Level: intermediate 2122 2123 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2124 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2125 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2126 2127 Notes: 2128 The number of rows and columns must be divisible by blocksize. 2129 This matrix type does not support complex Hermitian operation. 2130 2131 Specify the preallocated storage with either nz or nnz (not both). 2132 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2133 allocation. See Users-Manual: ch_mat for details. 2134 2135 If the nnz parameter is given then the nz parameter is ignored 2136 2137 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2138 @*/ 2139 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2140 { 2141 PetscErrorCode ierr; 2142 2143 PetscFunctionBegin; 2144 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2145 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2146 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2147 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2148 PetscFunctionReturn(0); 2149 } 2150 2151 #undef __FUNCT__ 2152 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2153 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2154 { 2155 Mat C; 2156 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2157 PetscErrorCode ierr; 2158 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2159 2160 PetscFunctionBegin; 2161 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2162 2163 *B = 0; 2164 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2165 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2166 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2167 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2168 c = (Mat_SeqSBAIJ*)C->data; 2169 2170 C->preallocated = PETSC_TRUE; 2171 C->factortype = A->factortype; 2172 c->row = 0; 2173 c->icol = 0; 2174 c->saved_values = 0; 2175 c->keepnonzeropattern = a->keepnonzeropattern; 2176 C->assembled = PETSC_TRUE; 2177 2178 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2179 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2180 c->bs2 = a->bs2; 2181 c->mbs = a->mbs; 2182 c->nbs = a->nbs; 2183 2184 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2185 c->imax = a->imax; 2186 c->ilen = a->ilen; 2187 c->free_imax_ilen = PETSC_FALSE; 2188 } else { 2189 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr); 2190 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2191 for (i=0; i<mbs; i++) { 2192 c->imax[i] = a->imax[i]; 2193 c->ilen[i] = a->ilen[i]; 2194 } 2195 c->free_imax_ilen = PETSC_TRUE; 2196 } 2197 2198 /* allocate the matrix space */ 2199 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2200 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr); 2201 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2202 c->i = a->i; 2203 c->j = a->j; 2204 c->singlemalloc = PETSC_FALSE; 2205 c->free_a = PETSC_TRUE; 2206 c->free_ij = PETSC_FALSE; 2207 c->parent = A; 2208 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2209 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2210 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2211 } else { 2212 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr); 2213 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2214 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2215 c->singlemalloc = PETSC_TRUE; 2216 c->free_a = PETSC_TRUE; 2217 c->free_ij = PETSC_TRUE; 2218 } 2219 if (mbs > 0) { 2220 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2221 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2222 } 2223 if (cpvalues == MAT_COPY_VALUES) { 2224 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2225 } else { 2226 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2227 } 2228 if (a->jshort) { 2229 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2230 /* if the parent matrix is reassembled, this child matrix will never notice */ 2231 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr); 2232 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2233 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2234 2235 c->free_jshort = PETSC_TRUE; 2236 } 2237 } 2238 2239 c->roworiented = a->roworiented; 2240 c->nonew = a->nonew; 2241 2242 if (a->diag) { 2243 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2244 c->diag = a->diag; 2245 c->free_diag = PETSC_FALSE; 2246 } else { 2247 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr); 2248 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2249 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2250 c->free_diag = PETSC_TRUE; 2251 } 2252 } 2253 c->nz = a->nz; 2254 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2255 c->solve_work = 0; 2256 c->mult_work = 0; 2257 2258 *B = C; 2259 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2260 PetscFunctionReturn(0); 2261 } 2262 2263 #undef __FUNCT__ 2264 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2265 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2266 { 2267 Mat_SeqSBAIJ *a; 2268 PetscErrorCode ierr; 2269 int fd; 2270 PetscMPIInt size; 2271 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs = newmat->rmap->bs; 2272 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2273 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2274 PetscInt *masked,nmask,tmp,bs2,ishift; 2275 PetscScalar *aa; 2276 MPI_Comm comm; 2277 2278 PetscFunctionBegin; 2279 /* force binary viewer to load .info file if it has not yet done so */ 2280 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 2281 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2282 ierr = PetscOptionsGetInt(((PetscObject)newmat)->options,((PetscObject)newmat)->prefix,"-matload_block_size",&bs,NULL);CHKERRQ(ierr); 2283 if (bs < 0) bs = 1; 2284 bs2 = bs*bs; 2285 2286 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2287 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2288 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2289 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2290 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2291 M = header[1]; N = header[2]; nz = header[3]; 2292 2293 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2294 2295 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2296 2297 /* 2298 This code adds extra rows to make sure the number of rows is 2299 divisible by the blocksize 2300 */ 2301 mbs = M/bs; 2302 extra_rows = bs - M + bs*(mbs); 2303 if (extra_rows == bs) extra_rows = 0; 2304 else mbs++; 2305 if (extra_rows) { 2306 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2307 } 2308 2309 /* Set global sizes if not already set */ 2310 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2311 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2312 } else { /* Check if the matrix global sizes are correct */ 2313 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2314 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); 2315 } 2316 2317 /* read in row lengths */ 2318 ierr = PetscMalloc1(M+extra_rows,&rowlengths);CHKERRQ(ierr); 2319 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2320 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2321 2322 /* read in column indices */ 2323 ierr = PetscMalloc1(nz+extra_rows,&jj);CHKERRQ(ierr); 2324 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2325 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2326 2327 /* loop over row lengths determining block row lengths */ 2328 ierr = PetscCalloc1(mbs,&s_browlengths);CHKERRQ(ierr); 2329 ierr = PetscMalloc2(mbs,&mask,mbs,&masked);CHKERRQ(ierr); 2330 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2331 rowcount = 0; 2332 nzcount = 0; 2333 for (i=0; i<mbs; i++) { 2334 nmask = 0; 2335 for (j=0; j<bs; j++) { 2336 kmax = rowlengths[rowcount]; 2337 for (k=0; k<kmax; k++) { 2338 tmp = jj[nzcount++]/bs; /* block col. index */ 2339 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2340 } 2341 rowcount++; 2342 } 2343 s_browlengths[i] += nmask; 2344 2345 /* zero out the mask elements we set */ 2346 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2347 } 2348 2349 /* Do preallocation */ 2350 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2351 a = (Mat_SeqSBAIJ*)newmat->data; 2352 2353 /* set matrix "i" values */ 2354 a->i[0] = 0; 2355 for (i=1; i<= mbs; i++) { 2356 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2357 a->ilen[i-1] = s_browlengths[i-1]; 2358 } 2359 a->nz = a->i[mbs]; 2360 2361 /* read in nonzero values */ 2362 ierr = PetscMalloc1(nz+extra_rows,&aa);CHKERRQ(ierr); 2363 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2364 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2365 2366 /* set "a" and "j" values into matrix */ 2367 nzcount = 0; jcount = 0; 2368 for (i=0; i<mbs; i++) { 2369 nzcountb = nzcount; 2370 nmask = 0; 2371 for (j=0; j<bs; j++) { 2372 kmax = rowlengths[i*bs+j]; 2373 for (k=0; k<kmax; k++) { 2374 tmp = jj[nzcount++]/bs; /* block col. index */ 2375 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2376 } 2377 } 2378 /* sort the masked values */ 2379 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2380 2381 /* set "j" values into matrix */ 2382 maskcount = 1; 2383 for (j=0; j<nmask; j++) { 2384 a->j[jcount++] = masked[j]; 2385 mask[masked[j]] = maskcount++; 2386 } 2387 2388 /* set "a" values into matrix */ 2389 ishift = bs2*a->i[i]; 2390 for (j=0; j<bs; j++) { 2391 kmax = rowlengths[i*bs+j]; 2392 for (k=0; k<kmax; k++) { 2393 tmp = jj[nzcountb]/bs; /* block col. index */ 2394 if (tmp >= i) { 2395 block = mask[tmp] - 1; 2396 point = jj[nzcountb] - bs*tmp; 2397 idx = ishift + bs2*block + j + bs*point; 2398 a->a[idx] = aa[nzcountb]; 2399 } 2400 nzcountb++; 2401 } 2402 } 2403 /* zero out the mask elements we set */ 2404 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2405 } 2406 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2407 2408 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2409 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2410 ierr = PetscFree(aa);CHKERRQ(ierr); 2411 ierr = PetscFree(jj);CHKERRQ(ierr); 2412 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2413 2414 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2415 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2416 PetscFunctionReturn(0); 2417 } 2418 2419 #undef __FUNCT__ 2420 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2421 /*@ 2422 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2423 (upper triangular entries in CSR format) provided by the user. 2424 2425 Collective on MPI_Comm 2426 2427 Input Parameters: 2428 + comm - must be an MPI communicator of size 1 2429 . bs - size of block 2430 . m - number of rows 2431 . n - number of columns 2432 . i - row indices 2433 . j - column indices 2434 - a - matrix values 2435 2436 Output Parameter: 2437 . mat - the matrix 2438 2439 Level: advanced 2440 2441 Notes: 2442 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2443 once the matrix is destroyed 2444 2445 You cannot set new nonzero locations into this matrix, that will generate an error. 2446 2447 The i and j indices are 0 based 2448 2449 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 2450 it is the regular CSR format excluding the lower triangular elements. 2451 2452 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2453 2454 @*/ 2455 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 2456 { 2457 PetscErrorCode ierr; 2458 PetscInt ii; 2459 Mat_SeqSBAIJ *sbaij; 2460 2461 PetscFunctionBegin; 2462 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2463 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2464 2465 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2466 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2467 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2468 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2469 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2470 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr); 2471 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2472 2473 sbaij->i = i; 2474 sbaij->j = j; 2475 sbaij->a = a; 2476 2477 sbaij->singlemalloc = PETSC_FALSE; 2478 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2479 sbaij->free_a = PETSC_FALSE; 2480 sbaij->free_ij = PETSC_FALSE; 2481 2482 for (ii=0; ii<m; ii++) { 2483 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2484 #if defined(PETSC_USE_DEBUG) 2485 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]); 2486 #endif 2487 } 2488 #if defined(PETSC_USE_DEBUG) 2489 for (ii=0; ii<sbaij->i[m]; ii++) { 2490 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2491 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]); 2492 } 2493 #endif 2494 2495 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2496 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2497 PetscFunctionReturn(0); 2498 } 2499 2500 #undef __FUNCT__ 2501 #define __FUNCT__ "MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ" 2502 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2503 { 2504 PetscErrorCode ierr; 2505 2506 PetscFunctionBegin; 2507 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2508 PetscFunctionReturn(0); 2509 } 2510 2511 2512 2513 2514