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