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