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