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