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