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