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