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