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