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