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