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