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