1 /* 2 Defines the basic matrix operations for the AIJ (compressed row) 3 matrix storage format. 4 */ 5 6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 #include <petscbt.h> 9 #include <petsc/private/kernels/blocktranspose.h> 10 11 PetscErrorCode MatSeqAIJSetTypeFromOptions(Mat A) 12 { 13 PetscBool flg; 14 char type[256]; 15 16 PetscFunctionBegin; 17 PetscObjectOptionsBegin((PetscObject)A); 18 PetscCall(PetscOptionsFList("-mat_seqaij_type","Matrix SeqAIJ type","MatSeqAIJSetType",MatSeqAIJList,"seqaij",type,256,&flg)); 19 if (flg) PetscCall(MatSeqAIJSetType(A,type)); 20 PetscOptionsEnd(); 21 PetscFunctionReturn(0); 22 } 23 24 PetscErrorCode MatGetColumnReductions_SeqAIJ(Mat A,PetscInt type,PetscReal *reductions) 25 { 26 PetscInt i,m,n; 27 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 28 29 PetscFunctionBegin; 30 PetscCall(MatGetSize(A,&m,&n)); 31 PetscCall(PetscArrayzero(reductions,n)); 32 if (type == NORM_2) { 33 for (i=0; i<aij->i[m]; i++) { 34 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 35 } 36 } else if (type == NORM_1) { 37 for (i=0; i<aij->i[m]; i++) { 38 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]); 39 } 40 } else if (type == NORM_INFINITY) { 41 for (i=0; i<aij->i[m]; i++) { 42 reductions[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),reductions[aij->j[i]]); 43 } 44 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 45 for (i=0; i<aij->i[m]; i++) { 46 reductions[aij->j[i]] += PetscRealPart(aij->a[i]); 47 } 48 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 49 for (i=0; i<aij->i[m]; i++) { 50 reductions[aij->j[i]] += PetscImaginaryPart(aij->a[i]); 51 } 52 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown reduction type"); 53 54 if (type == NORM_2) { 55 for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 56 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 57 for (i=0; i<n; i++) reductions[i] /= m; 58 } 59 PetscFunctionReturn(0); 60 } 61 62 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 63 { 64 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 65 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 66 const PetscInt *jj = a->j,*ii = a->i; 67 PetscInt *rows; 68 69 PetscFunctionBegin; 70 for (i=0; i<m; i++) { 71 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 72 cnt++; 73 } 74 } 75 PetscCall(PetscMalloc1(cnt,&rows)); 76 cnt = 0; 77 for (i=0; i<m; i++) { 78 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 79 rows[cnt] = i; 80 cnt++; 81 } 82 } 83 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is)); 84 PetscFunctionReturn(0); 85 } 86 87 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 88 { 89 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 90 const MatScalar *aa; 91 PetscInt i,m=A->rmap->n,cnt = 0; 92 const PetscInt *ii = a->i,*jj = a->j,*diag; 93 PetscInt *rows; 94 95 PetscFunctionBegin; 96 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 97 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 98 diag = a->diag; 99 for (i=0; i<m; i++) { 100 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 101 cnt++; 102 } 103 } 104 PetscCall(PetscMalloc1(cnt,&rows)); 105 cnt = 0; 106 for (i=0; i<m; i++) { 107 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 108 rows[cnt++] = i; 109 } 110 } 111 *nrows = cnt; 112 *zrows = rows; 113 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 118 { 119 PetscInt nrows,*rows; 120 121 PetscFunctionBegin; 122 *zrows = NULL; 123 PetscCall(MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows)); 124 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 129 { 130 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 131 const MatScalar *aa; 132 PetscInt m=A->rmap->n,cnt = 0; 133 const PetscInt *ii; 134 PetscInt n,i,j,*rows; 135 136 PetscFunctionBegin; 137 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 138 *keptrows = NULL; 139 ii = a->i; 140 for (i=0; i<m; i++) { 141 n = ii[i+1] - ii[i]; 142 if (!n) { 143 cnt++; 144 goto ok1; 145 } 146 for (j=ii[i]; j<ii[i+1]; j++) { 147 if (aa[j] != 0.0) goto ok1; 148 } 149 cnt++; 150 ok1:; 151 } 152 if (!cnt) { 153 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 154 PetscFunctionReturn(0); 155 } 156 PetscCall(PetscMalloc1(A->rmap->n-cnt,&rows)); 157 cnt = 0; 158 for (i=0; i<m; i++) { 159 n = ii[i+1] - ii[i]; 160 if (!n) continue; 161 for (j=ii[i]; j<ii[i+1]; j++) { 162 if (aa[j] != 0.0) { 163 rows[cnt++] = i; 164 break; 165 } 166 } 167 } 168 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 169 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows)); 170 PetscFunctionReturn(0); 171 } 172 173 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 174 { 175 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 176 PetscInt i,m = Y->rmap->n; 177 const PetscInt *diag; 178 MatScalar *aa; 179 const PetscScalar *v; 180 PetscBool missing; 181 182 PetscFunctionBegin; 183 if (Y->assembled) { 184 PetscCall(MatMissingDiagonal_SeqAIJ(Y,&missing,NULL)); 185 if (!missing) { 186 diag = aij->diag; 187 PetscCall(VecGetArrayRead(D,&v)); 188 PetscCall(MatSeqAIJGetArray(Y,&aa)); 189 if (is == INSERT_VALUES) { 190 for (i=0; i<m; i++) { 191 aa[diag[i]] = v[i]; 192 } 193 } else { 194 for (i=0; i<m; i++) { 195 aa[diag[i]] += v[i]; 196 } 197 } 198 PetscCall(MatSeqAIJRestoreArray(Y,&aa)); 199 PetscCall(VecRestoreArrayRead(D,&v)); 200 PetscFunctionReturn(0); 201 } 202 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 203 } 204 PetscCall(MatDiagonalSet_Default(Y,D,is)); 205 PetscFunctionReturn(0); 206 } 207 208 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 209 { 210 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 211 PetscInt i,ishift; 212 213 PetscFunctionBegin; 214 if (m) *m = A->rmap->n; 215 if (!ia) PetscFunctionReturn(0); 216 ishift = 0; 217 if (symmetric && A->structurally_symmetric != PETSC_BOOL3_TRUE) { 218 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja)); 219 } else if (oshift == 1) { 220 PetscInt *tia; 221 PetscInt nz = a->i[A->rmap->n]; 222 /* malloc space and add 1 to i and j indices */ 223 PetscCall(PetscMalloc1(A->rmap->n+1,&tia)); 224 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 225 *ia = tia; 226 if (ja) { 227 PetscInt *tja; 228 PetscCall(PetscMalloc1(nz+1,&tja)); 229 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 230 *ja = tja; 231 } 232 } else { 233 *ia = a->i; 234 if (ja) *ja = a->j; 235 } 236 PetscFunctionReturn(0); 237 } 238 239 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 240 { 241 PetscFunctionBegin; 242 if (!ia) PetscFunctionReturn(0); 243 if ((symmetric && A->structurally_symmetric != PETSC_BOOL3_TRUE) || oshift == 1) { 244 PetscCall(PetscFree(*ia)); 245 if (ja) PetscCall(PetscFree(*ja)); 246 } 247 PetscFunctionReturn(0); 248 } 249 250 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 251 { 252 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 253 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 254 PetscInt nz = a->i[m],row,*jj,mr,col; 255 256 PetscFunctionBegin; 257 *nn = n; 258 if (!ia) PetscFunctionReturn(0); 259 if (symmetric) { 260 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja)); 261 } else { 262 PetscCall(PetscCalloc1(n,&collengths)); 263 PetscCall(PetscMalloc1(n+1,&cia)); 264 PetscCall(PetscMalloc1(nz,&cja)); 265 jj = a->j; 266 for (i=0; i<nz; i++) { 267 collengths[jj[i]]++; 268 } 269 cia[0] = oshift; 270 for (i=0; i<n; i++) { 271 cia[i+1] = cia[i] + collengths[i]; 272 } 273 PetscCall(PetscArrayzero(collengths,n)); 274 jj = a->j; 275 for (row=0; row<m; row++) { 276 mr = a->i[row+1] - a->i[row]; 277 for (i=0; i<mr; i++) { 278 col = *jj++; 279 280 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 281 } 282 } 283 PetscCall(PetscFree(collengths)); 284 *ia = cia; *ja = cja; 285 } 286 PetscFunctionReturn(0); 287 } 288 289 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 290 { 291 PetscFunctionBegin; 292 if (!ia) PetscFunctionReturn(0); 293 294 PetscCall(PetscFree(*ia)); 295 PetscCall(PetscFree(*ja)); 296 PetscFunctionReturn(0); 297 } 298 299 /* 300 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 301 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 302 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 303 */ 304 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 305 { 306 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 307 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 308 PetscInt nz = a->i[m],row,mr,col,tmp; 309 PetscInt *cspidx; 310 const PetscInt *jj; 311 312 PetscFunctionBegin; 313 *nn = n; 314 if (!ia) PetscFunctionReturn(0); 315 316 PetscCall(PetscCalloc1(n,&collengths)); 317 PetscCall(PetscMalloc1(n+1,&cia)); 318 PetscCall(PetscMalloc1(nz,&cja)); 319 PetscCall(PetscMalloc1(nz,&cspidx)); 320 jj = a->j; 321 for (i=0; i<nz; i++) { 322 collengths[jj[i]]++; 323 } 324 cia[0] = oshift; 325 for (i=0; i<n; i++) { 326 cia[i+1] = cia[i] + collengths[i]; 327 } 328 PetscCall(PetscArrayzero(collengths,n)); 329 jj = a->j; 330 for (row=0; row<m; row++) { 331 mr = a->i[row+1] - a->i[row]; 332 for (i=0; i<mr; i++) { 333 col = *jj++; 334 tmp = cia[col] + collengths[col]++ - oshift; 335 cspidx[tmp] = a->i[row] + i; /* index of a->j */ 336 cja[tmp] = row + oshift; 337 } 338 } 339 PetscCall(PetscFree(collengths)); 340 *ia = cia; 341 *ja = cja; 342 *spidx = cspidx; 343 PetscFunctionReturn(0); 344 } 345 346 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 347 { 348 PetscFunctionBegin; 349 PetscCall(MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done)); 350 PetscCall(PetscFree(*spidx)); 351 PetscFunctionReturn(0); 352 } 353 354 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 355 { 356 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 357 PetscInt *ai = a->i; 358 PetscScalar *aa; 359 360 PetscFunctionBegin; 361 PetscCall(MatSeqAIJGetArray(A,&aa)); 362 PetscCall(PetscArraycpy(aa+ai[row],v,ai[row+1]-ai[row])); 363 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 364 PetscFunctionReturn(0); 365 } 366 367 /* 368 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 369 370 - a single row of values is set with each call 371 - no row or column indices are negative or (in error) larger than the number of rows or columns 372 - the values are always added to the matrix, not set 373 - no new locations are introduced in the nonzero structure of the matrix 374 375 This does NOT assume the global column indices are sorted 376 377 */ 378 379 #include <petsc/private/isimpl.h> 380 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 381 { 382 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 383 PetscInt low,high,t,row,nrow,i,col,l; 384 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 385 PetscInt lastcol = -1; 386 MatScalar *ap,value,*aa; 387 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 388 389 PetscFunctionBegin; 390 PetscCall(MatSeqAIJGetArray(A,&aa)); 391 row = ridx[im[0]]; 392 rp = aj + ai[row]; 393 ap = aa + ai[row]; 394 nrow = ailen[row]; 395 low = 0; 396 high = nrow; 397 for (l=0; l<n; l++) { /* loop over added columns */ 398 col = cidx[in[l]]; 399 value = v[l]; 400 401 if (col <= lastcol) low = 0; 402 else high = nrow; 403 lastcol = col; 404 while (high-low > 5) { 405 t = (low+high)/2; 406 if (rp[t] > col) high = t; 407 else low = t; 408 } 409 for (i=low; i<high; i++) { 410 if (rp[i] == col) { 411 ap[i] += value; 412 low = i + 1; 413 break; 414 } 415 } 416 } 417 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 418 return 0; 419 } 420 421 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 422 { 423 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 424 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 425 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 426 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 427 MatScalar *ap=NULL,value=0.0,*aa; 428 PetscBool ignorezeroentries = a->ignorezeroentries; 429 PetscBool roworiented = a->roworiented; 430 431 PetscFunctionBegin; 432 PetscCall(MatSeqAIJGetArray(A,&aa)); 433 for (k=0; k<m; k++) { /* loop over added rows */ 434 row = im[k]; 435 if (row < 0) continue; 436 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 437 rp = aj + ai[row]; 438 if (!A->structure_only) ap = aa + ai[row]; 439 rmax = imax[row]; nrow = ailen[row]; 440 low = 0; 441 high = nrow; 442 for (l=0; l<n; l++) { /* loop over added columns */ 443 if (in[l] < 0) continue; 444 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 445 col = in[l]; 446 if (v && !A->structure_only) value = roworiented ? v[l + k*n] : v[k + l*m]; 447 if (!A->structure_only && value == 0.0 && ignorezeroentries && is == ADD_VALUES && row != col) continue; 448 449 if (col <= lastcol) low = 0; 450 else high = nrow; 451 lastcol = col; 452 while (high-low > 5) { 453 t = (low+high)/2; 454 if (rp[t] > col) high = t; 455 else low = t; 456 } 457 for (i=low; i<high; i++) { 458 if (rp[i] > col) break; 459 if (rp[i] == col) { 460 if (!A->structure_only) { 461 if (is == ADD_VALUES) { 462 ap[i] += value; 463 (void)PetscLogFlops(1.0); 464 } 465 else ap[i] = value; 466 } 467 low = i + 1; 468 goto noinsert; 469 } 470 } 471 if (value == 0.0 && ignorezeroentries && row != col) goto noinsert; 472 if (nonew == 1) goto noinsert; 473 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%" PetscInt_FMT ",%" PetscInt_FMT ") in the matrix",row,col); 474 if (A->structure_only) { 475 MatSeqXAIJReallocateAIJ_structure_only(A,A->rmap->n,1,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar); 476 } else { 477 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 478 } 479 N = nrow++ - 1; a->nz++; high++; 480 /* shift up all the later entries in this row */ 481 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 482 rp[i] = col; 483 if (!A->structure_only) { 484 PetscCall(PetscArraymove(ap+i+1,ap+i,N-i+1)); 485 ap[i] = value; 486 } 487 low = i + 1; 488 A->nonzerostate++; 489 noinsert:; 490 } 491 ailen[row] = nrow; 492 } 493 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 494 PetscFunctionReturn(0); 495 } 496 497 PetscErrorCode MatSetValues_SeqAIJ_SortedFullNoPreallocation(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 498 { 499 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 500 PetscInt *rp,k,row; 501 PetscInt *ai = a->i; 502 PetscInt *aj = a->j; 503 MatScalar *aa,*ap; 504 505 PetscFunctionBegin; 506 PetscCheck(!A->was_assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot call on assembled matrix."); 507 PetscCheck(m*n+a->nz <= a->maxnz,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of entries in matrix will be larger than maximum nonzeros allocated for %" PetscInt_FMT " in MatSeqAIJSetTotalPreallocation()",a->maxnz); 508 509 PetscCall(MatSeqAIJGetArray(A,&aa)); 510 for (k=0; k<m; k++) { /* loop over added rows */ 511 row = im[k]; 512 rp = aj + ai[row]; 513 ap = aa + ai[row]; 514 515 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 516 if (!A->structure_only) { 517 if (v) { 518 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 519 v += n; 520 } else { 521 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 522 } 523 } 524 a->ilen[row] = n; 525 a->imax[row] = n; 526 a->i[row+1] = a->i[row]+n; 527 a->nz += n; 528 } 529 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 530 PetscFunctionReturn(0); 531 } 532 533 /*@ 534 MatSeqAIJSetTotalPreallocation - Sets an upper bound on the total number of expected nonzeros in the matrix. 535 536 Input Parameters: 537 + A - the SeqAIJ matrix 538 - nztotal - bound on the number of nonzeros 539 540 Level: advanced 541 542 Notes: 543 This can be called if you will be provided the matrix row by row (from row zero) with sorted column indices for each row. 544 Simply call MatSetValues() after this call to provide the matrix entries in the usual manner. This matrix may be used 545 as always with multiple matrix assemblies. 546 547 .seealso: `MatSetOption()`, `MAT_SORTED_FULL`, `MatSetValues()`, `MatSeqAIJSetPreallocation()` 548 @*/ 549 550 PetscErrorCode MatSeqAIJSetTotalPreallocation(Mat A,PetscInt nztotal) 551 { 552 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 553 554 PetscFunctionBegin; 555 PetscCall(PetscLayoutSetUp(A->rmap)); 556 PetscCall(PetscLayoutSetUp(A->cmap)); 557 a->maxnz = nztotal; 558 if (!a->imax) { 559 PetscCall(PetscMalloc1(A->rmap->n,&a->imax)); 560 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 561 } 562 if (!a->ilen) { 563 PetscCall(PetscMalloc1(A->rmap->n,&a->ilen)); 564 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 565 } else { 566 PetscCall(PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt))); 567 } 568 569 /* allocate the matrix space */ 570 if (A->structure_only) { 571 PetscCall(PetscMalloc1(nztotal,&a->j)); 572 PetscCall(PetscMalloc1(A->rmap->n+1,&a->i)); 573 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*sizeof(PetscInt))); 574 } else { 575 PetscCall(PetscMalloc3(nztotal,&a->a,nztotal,&a->j,A->rmap->n+1,&a->i)); 576 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*(sizeof(PetscScalar)+sizeof(PetscInt)))); 577 } 578 a->i[0] = 0; 579 if (A->structure_only) { 580 a->singlemalloc = PETSC_FALSE; 581 a->free_a = PETSC_FALSE; 582 } else { 583 a->singlemalloc = PETSC_TRUE; 584 a->free_a = PETSC_TRUE; 585 } 586 a->free_ij = PETSC_TRUE; 587 A->ops->setvalues = MatSetValues_SeqAIJ_SortedFullNoPreallocation; 588 A->preallocated = PETSC_TRUE; 589 PetscFunctionReturn(0); 590 } 591 592 PetscErrorCode MatSetValues_SeqAIJ_SortedFull(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 593 { 594 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 595 PetscInt *rp,k,row; 596 PetscInt *ai = a->i,*ailen = a->ilen; 597 PetscInt *aj = a->j; 598 MatScalar *aa,*ap; 599 600 PetscFunctionBegin; 601 PetscCall(MatSeqAIJGetArray(A,&aa)); 602 for (k=0; k<m; k++) { /* loop over added rows */ 603 row = im[k]; 604 PetscCheck(n <= a->imax[row],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Preallocation for row %" PetscInt_FMT " does not match number of columns provided",n); 605 rp = aj + ai[row]; 606 ap = aa + ai[row]; 607 if (!A->was_assembled) { 608 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 609 } 610 if (!A->structure_only) { 611 if (v) { 612 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 613 v += n; 614 } else { 615 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 616 } 617 } 618 ailen[row] = n; 619 a->nz += n; 620 } 621 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 622 PetscFunctionReturn(0); 623 } 624 625 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 626 { 627 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 628 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 629 PetscInt *ai = a->i,*ailen = a->ilen; 630 MatScalar *ap,*aa; 631 632 PetscFunctionBegin; 633 PetscCall(MatSeqAIJGetArray(A,&aa)); 634 for (k=0; k<m; k++) { /* loop over rows */ 635 row = im[k]; 636 if (row < 0) {v += n; continue;} /* negative row */ 637 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 638 rp = aj + ai[row]; ap = aa + ai[row]; 639 nrow = ailen[row]; 640 for (l=0; l<n; l++) { /* loop over columns */ 641 if (in[l] < 0) {v++; continue;} /* negative column */ 642 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 643 col = in[l]; 644 high = nrow; low = 0; /* assume unsorted */ 645 while (high-low > 5) { 646 t = (low+high)/2; 647 if (rp[t] > col) high = t; 648 else low = t; 649 } 650 for (i=low; i<high; i++) { 651 if (rp[i] > col) break; 652 if (rp[i] == col) { 653 *v++ = ap[i]; 654 goto finished; 655 } 656 } 657 *v++ = 0.0; 658 finished:; 659 } 660 } 661 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 662 PetscFunctionReturn(0); 663 } 664 665 PetscErrorCode MatView_SeqAIJ_Binary(Mat mat,PetscViewer viewer) 666 { 667 Mat_SeqAIJ *A = (Mat_SeqAIJ*)mat->data; 668 const PetscScalar *av; 669 PetscInt header[4],M,N,m,nz,i; 670 PetscInt *rowlens; 671 672 PetscFunctionBegin; 673 PetscCall(PetscViewerSetUp(viewer)); 674 675 M = mat->rmap->N; 676 N = mat->cmap->N; 677 m = mat->rmap->n; 678 nz = A->nz; 679 680 /* write matrix header */ 681 header[0] = MAT_FILE_CLASSID; 682 header[1] = M; header[2] = N; header[3] = nz; 683 PetscCall(PetscViewerBinaryWrite(viewer,header,4,PETSC_INT)); 684 685 /* fill in and store row lengths */ 686 PetscCall(PetscMalloc1(m,&rowlens)); 687 for (i=0; i<m; i++) rowlens[i] = A->i[i+1] - A->i[i]; 688 PetscCall(PetscViewerBinaryWrite(viewer,rowlens,m,PETSC_INT)); 689 PetscCall(PetscFree(rowlens)); 690 /* store column indices */ 691 PetscCall(PetscViewerBinaryWrite(viewer,A->j,nz,PETSC_INT)); 692 /* store nonzero values */ 693 PetscCall(MatSeqAIJGetArrayRead(mat,&av)); 694 PetscCall(PetscViewerBinaryWrite(viewer,av,nz,PETSC_SCALAR)); 695 PetscCall(MatSeqAIJRestoreArrayRead(mat,&av)); 696 697 /* write block size option to the viewer's .info file */ 698 PetscCall(MatView_Binary_BlockSizes(mat,viewer)); 699 PetscFunctionReturn(0); 700 } 701 702 static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A,PetscViewer viewer) 703 { 704 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 705 PetscInt i,k,m=A->rmap->N; 706 707 PetscFunctionBegin; 708 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 709 for (i=0; i<m; i++) { 710 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 711 for (k=a->i[i]; k<a->i[i+1]; k++) { 712 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ") ",a->j[k])); 713 } 714 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 715 } 716 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 717 PetscFunctionReturn(0); 718 } 719 720 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 721 722 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 723 { 724 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 725 const PetscScalar *av; 726 PetscInt i,j,m = A->rmap->n; 727 const char *name; 728 PetscViewerFormat format; 729 730 PetscFunctionBegin; 731 if (A->structure_only) { 732 PetscCall(MatView_SeqAIJ_ASCII_structonly(A,viewer)); 733 PetscFunctionReturn(0); 734 } 735 736 PetscCall(PetscViewerGetFormat(viewer,&format)); 737 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 738 739 /* trigger copy to CPU if needed */ 740 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 741 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 742 if (format == PETSC_VIEWER_ASCII_MATLAB) { 743 PetscInt nofinalvalue = 0; 744 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 745 /* Need a dummy value to ensure the dimension of the matrix. */ 746 nofinalvalue = 1; 747 } 748 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 749 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",m,A->cmap->n)); 750 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %" PetscInt_FMT " \n",a->nz)); 751 #if defined(PETSC_USE_COMPLEX) 752 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",4);\n",a->nz+nofinalvalue)); 753 #else 754 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",a->nz+nofinalvalue)); 755 #endif 756 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n")); 757 758 for (i=0; i<m; i++) { 759 for (j=a->i[i]; j<a->i[i+1]; j++) { 760 #if defined(PETSC_USE_COMPLEX) 761 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",i+1,a->j[j]+1,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 762 #else 763 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,a->j[j]+1,(double)a->a[j])); 764 #endif 765 } 766 } 767 if (nofinalvalue) { 768 #if defined(PETSC_USE_COMPLEX) 769 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",m,A->cmap->n,0.,0.)); 770 #else 771 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",m,A->cmap->n,0.0)); 772 #endif 773 } 774 PetscCall(PetscObjectGetName((PetscObject)A,&name)); 775 PetscCall(PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name)); 776 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 777 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 778 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 779 for (i=0; i<m; i++) { 780 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 781 for (j=a->i[i]; j<a->i[i+1]; j++) { 782 #if defined(PETSC_USE_COMPLEX) 783 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 784 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 785 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 786 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 787 } else if (PetscRealPart(a->a[j]) != 0.0) { 788 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 789 } 790 #else 791 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 792 #endif 793 } 794 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 795 } 796 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 797 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 798 PetscInt nzd=0,fshift=1,*sptr; 799 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 800 PetscCall(PetscMalloc1(m+1,&sptr)); 801 for (i=0; i<m; i++) { 802 sptr[i] = nzd+1; 803 for (j=a->i[i]; j<a->i[i+1]; j++) { 804 if (a->j[j] >= i) { 805 #if defined(PETSC_USE_COMPLEX) 806 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 807 #else 808 if (a->a[j] != 0.0) nzd++; 809 #endif 810 } 811 } 812 } 813 sptr[m] = nzd+1; 814 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n\n",m,nzd)); 815 for (i=0; i<m+1; i+=6) { 816 if (i+4<m) { 817 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5])); 818 } else if (i+3<m) { 819 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4])); 820 } else if (i+2<m) { 821 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3])); 822 } else if (i+1<m) { 823 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2])); 824 } else if (i<m) { 825 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1])); 826 } else { 827 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT "\n",sptr[i])); 828 } 829 } 830 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 831 PetscCall(PetscFree(sptr)); 832 for (i=0; i<m; i++) { 833 for (j=a->i[i]; j<a->i[i+1]; j++) { 834 if (a->j[j] >= i) PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " ",a->j[j]+fshift)); 835 } 836 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 837 } 838 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 839 for (i=0; i<m; i++) { 840 for (j=a->i[i]; j<a->i[i+1]; j++) { 841 if (a->j[j] >= i) { 842 #if defined(PETSC_USE_COMPLEX) 843 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 844 PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 845 } 846 #else 847 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j])); 848 #endif 849 } 850 } 851 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 852 } 853 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 854 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 855 PetscInt cnt = 0,jcnt; 856 PetscScalar value; 857 #if defined(PETSC_USE_COMPLEX) 858 PetscBool realonly = PETSC_TRUE; 859 860 for (i=0; i<a->i[m]; i++) { 861 if (PetscImaginaryPart(a->a[i]) != 0.0) { 862 realonly = PETSC_FALSE; 863 break; 864 } 865 } 866 #endif 867 868 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 869 for (i=0; i<m; i++) { 870 jcnt = 0; 871 for (j=0; j<A->cmap->n; j++) { 872 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 873 value = a->a[cnt++]; 874 jcnt++; 875 } else { 876 value = 0.0; 877 } 878 #if defined(PETSC_USE_COMPLEX) 879 if (realonly) { 880 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value))); 881 } else { 882 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value))); 883 } 884 #else 885 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value)); 886 #endif 887 } 888 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 889 } 890 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 891 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 892 PetscInt fshift=1; 893 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 894 #if defined(PETSC_USE_COMPLEX) 895 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n")); 896 #else 897 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n")); 898 #endif 899 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, A->cmap->n, a->nz)); 900 for (i=0; i<m; i++) { 901 for (j=a->i[i]; j<a->i[i+1]; j++) { 902 #if defined(PETSC_USE_COMPLEX) 903 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 904 #else 905 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j])); 906 #endif 907 } 908 } 909 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 910 } else { 911 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 912 if (A->factortype) { 913 for (i=0; i<m; i++) { 914 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 915 /* L part */ 916 for (j=a->i[i]; j<a->i[i+1]; j++) { 917 #if defined(PETSC_USE_COMPLEX) 918 if (PetscImaginaryPart(a->a[j]) > 0.0) { 919 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 920 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 921 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 922 } else { 923 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 924 } 925 #else 926 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 927 #endif 928 } 929 /* diagonal */ 930 j = a->diag[i]; 931 #if defined(PETSC_USE_COMPLEX) 932 if (PetscImaginaryPart(a->a[j]) > 0.0) { 933 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)PetscImaginaryPart(1.0/a->a[j]))); 934 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 935 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)(-PetscImaginaryPart(1.0/a->a[j])))); 936 } else { 937 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]))); 938 } 939 #else 940 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)(1.0/a->a[j]))); 941 #endif 942 943 /* U part */ 944 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 945 #if defined(PETSC_USE_COMPLEX) 946 if (PetscImaginaryPart(a->a[j]) > 0.0) { 947 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 948 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 949 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 950 } else { 951 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 952 } 953 #else 954 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 955 #endif 956 } 957 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 958 } 959 } else { 960 for (i=0; i<m; i++) { 961 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 962 for (j=a->i[i]; j<a->i[i+1]; j++) { 963 #if defined(PETSC_USE_COMPLEX) 964 if (PetscImaginaryPart(a->a[j]) > 0.0) { 965 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 966 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 967 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 968 } else { 969 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 970 } 971 #else 972 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 973 #endif 974 } 975 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 976 } 977 } 978 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 979 } 980 PetscCall(PetscViewerFlush(viewer)); 981 PetscFunctionReturn(0); 982 } 983 984 #include <petscdraw.h> 985 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 986 { 987 Mat A = (Mat) Aa; 988 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 989 PetscInt i,j,m = A->rmap->n; 990 int color; 991 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 992 PetscViewer viewer; 993 PetscViewerFormat format; 994 const PetscScalar *aa; 995 996 PetscFunctionBegin; 997 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 998 PetscCall(PetscViewerGetFormat(viewer,&format)); 999 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 1000 1001 /* loop over matrix elements drawing boxes */ 1002 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1003 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1004 PetscDrawCollectiveBegin(draw); 1005 /* Blue for negative, Cyan for zero and Red for positive */ 1006 color = PETSC_DRAW_BLUE; 1007 for (i=0; i<m; i++) { 1008 y_l = m - i - 1.0; y_r = y_l + 1.0; 1009 for (j=a->i[i]; j<a->i[i+1]; j++) { 1010 x_l = a->j[j]; x_r = x_l + 1.0; 1011 if (PetscRealPart(aa[j]) >= 0.) continue; 1012 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1013 } 1014 } 1015 color = PETSC_DRAW_CYAN; 1016 for (i=0; i<m; i++) { 1017 y_l = m - i - 1.0; y_r = y_l + 1.0; 1018 for (j=a->i[i]; j<a->i[i+1]; j++) { 1019 x_l = a->j[j]; x_r = x_l + 1.0; 1020 if (aa[j] != 0.) continue; 1021 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1022 } 1023 } 1024 color = PETSC_DRAW_RED; 1025 for (i=0; i<m; i++) { 1026 y_l = m - i - 1.0; y_r = y_l + 1.0; 1027 for (j=a->i[i]; j<a->i[i+1]; j++) { 1028 x_l = a->j[j]; x_r = x_l + 1.0; 1029 if (PetscRealPart(aa[j]) <= 0.) continue; 1030 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1031 } 1032 } 1033 PetscDrawCollectiveEnd(draw); 1034 } else { 1035 /* use contour shading to indicate magnitude of values */ 1036 /* first determine max of all nonzero values */ 1037 PetscReal minv = 0.0, maxv = 0.0; 1038 PetscInt nz = a->nz, count = 0; 1039 PetscDraw popup; 1040 1041 for (i=0; i<nz; i++) { 1042 if (PetscAbsScalar(aa[i]) > maxv) maxv = PetscAbsScalar(aa[i]); 1043 } 1044 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1045 PetscCall(PetscDrawGetPopup(draw,&popup)); 1046 PetscCall(PetscDrawScalePopup(popup,minv,maxv)); 1047 1048 PetscDrawCollectiveBegin(draw); 1049 for (i=0; i<m; i++) { 1050 y_l = m - i - 1.0; 1051 y_r = y_l + 1.0; 1052 for (j=a->i[i]; j<a->i[i+1]; j++) { 1053 x_l = a->j[j]; 1054 x_r = x_l + 1.0; 1055 color = PetscDrawRealToColor(PetscAbsScalar(aa[count]),minv,maxv); 1056 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1057 count++; 1058 } 1059 } 1060 PetscDrawCollectiveEnd(draw); 1061 } 1062 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #include <petscdraw.h> 1067 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 1068 { 1069 PetscDraw draw; 1070 PetscReal xr,yr,xl,yl,h,w; 1071 PetscBool isnull; 1072 1073 PetscFunctionBegin; 1074 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 1075 PetscCall(PetscDrawIsNull(draw,&isnull)); 1076 if (isnull) PetscFunctionReturn(0); 1077 1078 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1079 xr += w; yr += h; xl = -w; yl = -h; 1080 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 1081 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 1082 PetscCall(PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A)); 1083 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 1084 PetscCall(PetscDrawSave(draw)); 1085 PetscFunctionReturn(0); 1086 } 1087 1088 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 1089 { 1090 PetscBool iascii,isbinary,isdraw; 1091 1092 PetscFunctionBegin; 1093 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1094 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 1095 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 1096 if (iascii) PetscCall(MatView_SeqAIJ_ASCII(A,viewer)); 1097 else if (isbinary) PetscCall(MatView_SeqAIJ_Binary(A,viewer)); 1098 else if (isdraw) PetscCall(MatView_SeqAIJ_Draw(A,viewer)); 1099 PetscCall(MatView_SeqAIJ_Inode(A,viewer)); 1100 PetscFunctionReturn(0); 1101 } 1102 1103 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 1104 { 1105 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1106 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 1107 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 1108 MatScalar *aa = a->a,*ap; 1109 PetscReal ratio = 0.6; 1110 1111 PetscFunctionBegin; 1112 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 1113 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1114 if (A->was_assembled && A->ass_nonzerostate == A->nonzerostate) { 1115 /* we need to respect users asking to use or not the inodes routine in between matrix assemblies */ 1116 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1117 PetscFunctionReturn(0); 1118 } 1119 1120 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 1121 for (i=1; i<m; i++) { 1122 /* move each row back by the amount of empty slots (fshift) before it*/ 1123 fshift += imax[i-1] - ailen[i-1]; 1124 rmax = PetscMax(rmax,ailen[i]); 1125 if (fshift) { 1126 ip = aj + ai[i]; 1127 ap = aa + ai[i]; 1128 N = ailen[i]; 1129 PetscCall(PetscArraymove(ip-fshift,ip,N)); 1130 if (!A->structure_only) { 1131 PetscCall(PetscArraymove(ap-fshift,ap,N)); 1132 } 1133 } 1134 ai[i] = ai[i-1] + ailen[i-1]; 1135 } 1136 if (m) { 1137 fshift += imax[m-1] - ailen[m-1]; 1138 ai[m] = ai[m-1] + ailen[m-1]; 1139 } 1140 1141 /* reset ilen and imax for each row */ 1142 a->nonzerorowcnt = 0; 1143 if (A->structure_only) { 1144 PetscCall(PetscFree(a->imax)); 1145 PetscCall(PetscFree(a->ilen)); 1146 } else { /* !A->structure_only */ 1147 for (i=0; i<m; i++) { 1148 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1149 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1150 } 1151 } 1152 a->nz = ai[m]; 1153 PetscCheck(!fshift || a->nounused != -1,PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, fshift); 1154 1155 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1156 PetscCall(PetscInfo(A,"Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded,%" PetscInt_FMT " used\n",m,A->cmap->n,fshift,a->nz)); 1157 PetscCall(PetscInfo(A,"Number of mallocs during MatSetValues() is %" PetscInt_FMT "\n",a->reallocs)); 1158 PetscCall(PetscInfo(A,"Maximum nonzeros in any row is %" PetscInt_FMT "\n",rmax)); 1159 1160 A->info.mallocs += a->reallocs; 1161 a->reallocs = 0; 1162 A->info.nz_unneeded = (PetscReal)fshift; 1163 a->rmax = rmax; 1164 1165 if (!A->structure_only) { 1166 PetscCall(MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio)); 1167 } 1168 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1169 PetscFunctionReturn(0); 1170 } 1171 1172 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1173 { 1174 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1175 PetscInt i,nz = a->nz; 1176 MatScalar *aa; 1177 1178 PetscFunctionBegin; 1179 PetscCall(MatSeqAIJGetArray(A,&aa)); 1180 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1181 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1182 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1183 PetscFunctionReturn(0); 1184 } 1185 1186 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1187 { 1188 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1189 PetscInt i,nz = a->nz; 1190 MatScalar *aa; 1191 1192 PetscFunctionBegin; 1193 PetscCall(MatSeqAIJGetArray(A,&aa)); 1194 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1195 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1196 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1197 PetscFunctionReturn(0); 1198 } 1199 1200 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1201 { 1202 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1203 MatScalar *aa; 1204 1205 PetscFunctionBegin; 1206 PetscCall(MatSeqAIJGetArrayWrite(A,&aa)); 1207 PetscCall(PetscArrayzero(aa,a->i[A->rmap->n])); 1208 PetscCall(MatSeqAIJRestoreArrayWrite(A,&aa)); 1209 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1210 PetscFunctionReturn(0); 1211 } 1212 1213 PETSC_INTERN PetscErrorCode MatResetPreallocationCOO_SeqAIJ(Mat A) 1214 { 1215 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1216 1217 PetscFunctionBegin; 1218 PetscCall(PetscFree(a->perm)); 1219 PetscCall(PetscFree(a->jmap)); 1220 PetscFunctionReturn(0); 1221 } 1222 1223 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1224 { 1225 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1226 1227 PetscFunctionBegin; 1228 #if defined(PETSC_USE_LOG) 1229 PetscLogObjectState((PetscObject)A,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT,A->rmap->n,A->cmap->n,a->nz); 1230 #endif 1231 PetscCall(MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i)); 1232 PetscCall(MatResetPreallocationCOO_SeqAIJ(A)); 1233 PetscCall(ISDestroy(&a->row)); 1234 PetscCall(ISDestroy(&a->col)); 1235 PetscCall(PetscFree(a->diag)); 1236 PetscCall(PetscFree(a->ibdiag)); 1237 PetscCall(PetscFree(a->imax)); 1238 PetscCall(PetscFree(a->ilen)); 1239 PetscCall(PetscFree(a->ipre)); 1240 PetscCall(PetscFree3(a->idiag,a->mdiag,a->ssor_work)); 1241 PetscCall(PetscFree(a->solve_work)); 1242 PetscCall(ISDestroy(&a->icol)); 1243 PetscCall(PetscFree(a->saved_values)); 1244 PetscCall(PetscFree2(a->compressedrow.i,a->compressedrow.rindex)); 1245 PetscCall(MatDestroy_SeqAIJ_Inode(A)); 1246 PetscCall(PetscFree(A->data)); 1247 1248 /* MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted may allocate this. 1249 That function is so heavily used (sometimes in an hidden way through multnumeric function pointers) 1250 that is hard to properly add this data to the MatProduct data. We free it here to avoid 1251 users reusing the matrix object with different data to incur in obscure segmentation faults 1252 due to different matrix sizes */ 1253 PetscCall(PetscObjectCompose((PetscObject)A,"__PETSc__ab_dense",NULL)); 1254 1255 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 1256 PetscCall(PetscObjectComposeFunction((PetscObject)A,"PetscMatlabEnginePut_C",NULL)); 1257 PetscCall(PetscObjectComposeFunction((PetscObject)A,"PetscMatlabEngineGet_C",NULL)); 1258 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL)); 1259 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL)); 1260 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL)); 1261 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL)); 1262 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL)); 1263 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL)); 1264 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijsell_C",NULL)); 1265 #if defined(PETSC_HAVE_MKL_SPARSE) 1266 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijmkl_C",NULL)); 1267 #endif 1268 #if defined(PETSC_HAVE_CUDA) 1269 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcusparse_C",NULL)); 1270 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",NULL)); 1271 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",NULL)); 1272 #endif 1273 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 1274 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijkokkos_C",NULL)); 1275 #endif 1276 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcrl_C",NULL)); 1277 #if defined(PETSC_HAVE_ELEMENTAL) 1278 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL)); 1279 #endif 1280 #if defined(PETSC_HAVE_SCALAPACK) 1281 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_scalapack_C",NULL)); 1282 #endif 1283 #if defined(PETSC_HAVE_HYPRE) 1284 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL)); 1285 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",NULL)); 1286 #endif 1287 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL)); 1288 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsell_C",NULL)); 1289 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_is_C",NULL)); 1290 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL)); 1291 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatIsHermitianTranspose_C",NULL)); 1292 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL)); 1293 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL)); 1294 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL)); 1295 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL)); 1296 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_is_seqaij_C",NULL)); 1297 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqdense_seqaij_C",NULL)); 1298 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaij_C",NULL)); 1299 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJKron_C",NULL)); 1300 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL)); 1301 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL)); 1302 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL)); 1303 /* these calls do not belong here: the subclasses Duplicate/Destroy are wrong */ 1304 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaijsell_seqaij_C",NULL)); 1305 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaijperm_seqaij_C",NULL)); 1306 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijviennacl_C",NULL)); 1307 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijviennacl_seqdense_C",NULL)); 1308 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijviennacl_seqaij_C",NULL)); 1309 PetscFunctionReturn(0); 1310 } 1311 1312 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1313 { 1314 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1315 1316 PetscFunctionBegin; 1317 switch (op) { 1318 case MAT_ROW_ORIENTED: 1319 a->roworiented = flg; 1320 break; 1321 case MAT_KEEP_NONZERO_PATTERN: 1322 a->keepnonzeropattern = flg; 1323 break; 1324 case MAT_NEW_NONZERO_LOCATIONS: 1325 a->nonew = (flg ? 0 : 1); 1326 break; 1327 case MAT_NEW_NONZERO_LOCATION_ERR: 1328 a->nonew = (flg ? -1 : 0); 1329 break; 1330 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1331 a->nonew = (flg ? -2 : 0); 1332 break; 1333 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1334 a->nounused = (flg ? -1 : 0); 1335 break; 1336 case MAT_IGNORE_ZERO_ENTRIES: 1337 a->ignorezeroentries = flg; 1338 break; 1339 case MAT_SPD: 1340 case MAT_SYMMETRIC: 1341 case MAT_STRUCTURALLY_SYMMETRIC: 1342 case MAT_HERMITIAN: 1343 case MAT_SYMMETRY_ETERNAL: 1344 case MAT_STRUCTURE_ONLY: 1345 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 1346 case MAT_SPD_ETERNAL: 1347 /* if the diagonal matrix is square it inherits some of the properties above */ 1348 break; 1349 case MAT_FORCE_DIAGONAL_ENTRIES: 1350 case MAT_IGNORE_OFF_PROC_ENTRIES: 1351 case MAT_USE_HASH_TABLE: 1352 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 1353 break; 1354 case MAT_USE_INODES: 1355 PetscCall(MatSetOption_SeqAIJ_Inode(A,MAT_USE_INODES,flg)); 1356 break; 1357 case MAT_SUBMAT_SINGLEIS: 1358 A->submat_singleis = flg; 1359 break; 1360 case MAT_SORTED_FULL: 1361 if (flg) A->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 1362 else A->ops->setvalues = MatSetValues_SeqAIJ; 1363 break; 1364 case MAT_FORM_EXPLICIT_TRANSPOSE: 1365 A->form_explicit_transpose = flg; 1366 break; 1367 default: 1368 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1369 } 1370 PetscFunctionReturn(0); 1371 } 1372 1373 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1374 { 1375 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1376 PetscInt i,j,n,*ai=a->i,*aj=a->j; 1377 PetscScalar *x; 1378 const PetscScalar *aa; 1379 1380 PetscFunctionBegin; 1381 PetscCall(VecGetLocalSize(v,&n)); 1382 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1383 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1384 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1385 PetscInt *diag=a->diag; 1386 PetscCall(VecGetArrayWrite(v,&x)); 1387 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1388 PetscCall(VecRestoreArrayWrite(v,&x)); 1389 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1390 PetscFunctionReturn(0); 1391 } 1392 1393 PetscCall(VecGetArrayWrite(v,&x)); 1394 for (i=0; i<n; i++) { 1395 x[i] = 0.0; 1396 for (j=ai[i]; j<ai[i+1]; j++) { 1397 if (aj[j] == i) { 1398 x[i] = aa[j]; 1399 break; 1400 } 1401 } 1402 } 1403 PetscCall(VecRestoreArrayWrite(v,&x)); 1404 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1405 PetscFunctionReturn(0); 1406 } 1407 1408 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1409 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1410 { 1411 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1412 const MatScalar *aa; 1413 PetscScalar *y; 1414 const PetscScalar *x; 1415 PetscInt m = A->rmap->n; 1416 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1417 const MatScalar *v; 1418 PetscScalar alpha; 1419 PetscInt n,i,j; 1420 const PetscInt *idx,*ii,*ridx=NULL; 1421 Mat_CompressedRow cprow = a->compressedrow; 1422 PetscBool usecprow = cprow.use; 1423 #endif 1424 1425 PetscFunctionBegin; 1426 if (zz != yy) PetscCall(VecCopy(zz,yy)); 1427 PetscCall(VecGetArrayRead(xx,&x)); 1428 PetscCall(VecGetArray(yy,&y)); 1429 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1430 1431 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1432 fortranmulttransposeaddaij_(&m,x,a->i,a->j,aa,y); 1433 #else 1434 if (usecprow) { 1435 m = cprow.nrows; 1436 ii = cprow.i; 1437 ridx = cprow.rindex; 1438 } else { 1439 ii = a->i; 1440 } 1441 for (i=0; i<m; i++) { 1442 idx = a->j + ii[i]; 1443 v = aa + ii[i]; 1444 n = ii[i+1] - ii[i]; 1445 if (usecprow) { 1446 alpha = x[ridx[i]]; 1447 } else { 1448 alpha = x[i]; 1449 } 1450 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1451 } 1452 #endif 1453 PetscCall(PetscLogFlops(2.0*a->nz)); 1454 PetscCall(VecRestoreArrayRead(xx,&x)); 1455 PetscCall(VecRestoreArray(yy,&y)); 1456 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1457 PetscFunctionReturn(0); 1458 } 1459 1460 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1461 { 1462 PetscFunctionBegin; 1463 PetscCall(VecSet(yy,0.0)); 1464 PetscCall(MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy)); 1465 PetscFunctionReturn(0); 1466 } 1467 1468 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1469 1470 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1471 { 1472 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1473 PetscScalar *y; 1474 const PetscScalar *x; 1475 const MatScalar *aa,*a_a; 1476 PetscInt m=A->rmap->n; 1477 const PetscInt *aj,*ii,*ridx=NULL; 1478 PetscInt n,i; 1479 PetscScalar sum; 1480 PetscBool usecprow=a->compressedrow.use; 1481 1482 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1483 #pragma disjoint(*x,*y,*aa) 1484 #endif 1485 1486 PetscFunctionBegin; 1487 if (a->inode.use && a->inode.checked) { 1488 PetscCall(MatMult_SeqAIJ_Inode(A,xx,yy)); 1489 PetscFunctionReturn(0); 1490 } 1491 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1492 PetscCall(VecGetArrayRead(xx,&x)); 1493 PetscCall(VecGetArray(yy,&y)); 1494 ii = a->i; 1495 if (usecprow) { /* use compressed row format */ 1496 PetscCall(PetscArrayzero(y,m)); 1497 m = a->compressedrow.nrows; 1498 ii = a->compressedrow.i; 1499 ridx = a->compressedrow.rindex; 1500 for (i=0; i<m; i++) { 1501 n = ii[i+1] - ii[i]; 1502 aj = a->j + ii[i]; 1503 aa = a_a + ii[i]; 1504 sum = 0.0; 1505 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1506 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1507 y[*ridx++] = sum; 1508 } 1509 } else { /* do not use compressed row format */ 1510 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1511 aj = a->j; 1512 aa = a_a; 1513 fortranmultaij_(&m,x,ii,aj,aa,y); 1514 #else 1515 for (i=0; i<m; i++) { 1516 n = ii[i+1] - ii[i]; 1517 aj = a->j + ii[i]; 1518 aa = a_a + ii[i]; 1519 sum = 0.0; 1520 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1521 y[i] = sum; 1522 } 1523 #endif 1524 } 1525 PetscCall(PetscLogFlops(2.0*a->nz - a->nonzerorowcnt)); 1526 PetscCall(VecRestoreArrayRead(xx,&x)); 1527 PetscCall(VecRestoreArray(yy,&y)); 1528 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1529 PetscFunctionReturn(0); 1530 } 1531 1532 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1533 { 1534 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1535 PetscScalar *y; 1536 const PetscScalar *x; 1537 const MatScalar *aa,*a_a; 1538 PetscInt m=A->rmap->n; 1539 const PetscInt *aj,*ii,*ridx=NULL; 1540 PetscInt n,i,nonzerorow=0; 1541 PetscScalar sum; 1542 PetscBool usecprow=a->compressedrow.use; 1543 1544 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1545 #pragma disjoint(*x,*y,*aa) 1546 #endif 1547 1548 PetscFunctionBegin; 1549 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1550 PetscCall(VecGetArrayRead(xx,&x)); 1551 PetscCall(VecGetArray(yy,&y)); 1552 if (usecprow) { /* use compressed row format */ 1553 m = a->compressedrow.nrows; 1554 ii = a->compressedrow.i; 1555 ridx = a->compressedrow.rindex; 1556 for (i=0; i<m; i++) { 1557 n = ii[i+1] - ii[i]; 1558 aj = a->j + ii[i]; 1559 aa = a_a + ii[i]; 1560 sum = 0.0; 1561 nonzerorow += (n>0); 1562 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1563 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1564 y[*ridx++] = sum; 1565 } 1566 } else { /* do not use compressed row format */ 1567 ii = a->i; 1568 for (i=0; i<m; i++) { 1569 n = ii[i+1] - ii[i]; 1570 aj = a->j + ii[i]; 1571 aa = a_a + ii[i]; 1572 sum = 0.0; 1573 nonzerorow += (n>0); 1574 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1575 y[i] = sum; 1576 } 1577 } 1578 PetscCall(PetscLogFlops(2.0*a->nz - nonzerorow)); 1579 PetscCall(VecRestoreArrayRead(xx,&x)); 1580 PetscCall(VecRestoreArray(yy,&y)); 1581 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1582 PetscFunctionReturn(0); 1583 } 1584 1585 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1586 { 1587 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1588 PetscScalar *y,*z; 1589 const PetscScalar *x; 1590 const MatScalar *aa,*a_a; 1591 PetscInt m = A->rmap->n,*aj,*ii; 1592 PetscInt n,i,*ridx=NULL; 1593 PetscScalar sum; 1594 PetscBool usecprow=a->compressedrow.use; 1595 1596 PetscFunctionBegin; 1597 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1598 PetscCall(VecGetArrayRead(xx,&x)); 1599 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1600 if (usecprow) { /* use compressed row format */ 1601 if (zz != yy) { 1602 PetscCall(PetscArraycpy(z,y,m)); 1603 } 1604 m = a->compressedrow.nrows; 1605 ii = a->compressedrow.i; 1606 ridx = a->compressedrow.rindex; 1607 for (i=0; i<m; i++) { 1608 n = ii[i+1] - ii[i]; 1609 aj = a->j + ii[i]; 1610 aa = a_a + ii[i]; 1611 sum = y[*ridx]; 1612 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1613 z[*ridx++] = sum; 1614 } 1615 } else { /* do not use compressed row format */ 1616 ii = a->i; 1617 for (i=0; i<m; i++) { 1618 n = ii[i+1] - ii[i]; 1619 aj = a->j + ii[i]; 1620 aa = a_a + ii[i]; 1621 sum = y[i]; 1622 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1623 z[i] = sum; 1624 } 1625 } 1626 PetscCall(PetscLogFlops(2.0*a->nz)); 1627 PetscCall(VecRestoreArrayRead(xx,&x)); 1628 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1629 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1630 PetscFunctionReturn(0); 1631 } 1632 1633 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1634 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1635 { 1636 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1637 PetscScalar *y,*z; 1638 const PetscScalar *x; 1639 const MatScalar *aa,*a_a; 1640 const PetscInt *aj,*ii,*ridx=NULL; 1641 PetscInt m = A->rmap->n,n,i; 1642 PetscScalar sum; 1643 PetscBool usecprow=a->compressedrow.use; 1644 1645 PetscFunctionBegin; 1646 if (a->inode.use && a->inode.checked) { 1647 PetscCall(MatMultAdd_SeqAIJ_Inode(A,xx,yy,zz)); 1648 PetscFunctionReturn(0); 1649 } 1650 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1651 PetscCall(VecGetArrayRead(xx,&x)); 1652 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1653 if (usecprow) { /* use compressed row format */ 1654 if (zz != yy) { 1655 PetscCall(PetscArraycpy(z,y,m)); 1656 } 1657 m = a->compressedrow.nrows; 1658 ii = a->compressedrow.i; 1659 ridx = a->compressedrow.rindex; 1660 for (i=0; i<m; i++) { 1661 n = ii[i+1] - ii[i]; 1662 aj = a->j + ii[i]; 1663 aa = a_a + ii[i]; 1664 sum = y[*ridx]; 1665 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1666 z[*ridx++] = sum; 1667 } 1668 } else { /* do not use compressed row format */ 1669 ii = a->i; 1670 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1671 aj = a->j; 1672 aa = a_a; 1673 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1674 #else 1675 for (i=0; i<m; i++) { 1676 n = ii[i+1] - ii[i]; 1677 aj = a->j + ii[i]; 1678 aa = a_a + ii[i]; 1679 sum = y[i]; 1680 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1681 z[i] = sum; 1682 } 1683 #endif 1684 } 1685 PetscCall(PetscLogFlops(2.0*a->nz)); 1686 PetscCall(VecRestoreArrayRead(xx,&x)); 1687 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1688 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1689 PetscFunctionReturn(0); 1690 } 1691 1692 /* 1693 Adds diagonal pointers to sparse matrix structure. 1694 */ 1695 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1696 { 1697 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1698 PetscInt i,j,m = A->rmap->n; 1699 PetscBool alreadySet = PETSC_TRUE; 1700 1701 PetscFunctionBegin; 1702 if (!a->diag) { 1703 PetscCall(PetscMalloc1(m,&a->diag)); 1704 PetscCall(PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt))); 1705 alreadySet = PETSC_FALSE; 1706 } 1707 for (i=0; i<A->rmap->n; i++) { 1708 /* If A's diagonal is already correctly set, this fast track enables cheap and repeated MatMarkDiagonal_SeqAIJ() calls */ 1709 if (alreadySet) { 1710 PetscInt pos = a->diag[i]; 1711 if (pos >= a->i[i] && pos < a->i[i+1] && a->j[pos] == i) continue; 1712 } 1713 1714 a->diag[i] = a->i[i+1]; 1715 for (j=a->i[i]; j<a->i[i+1]; j++) { 1716 if (a->j[j] == i) { 1717 a->diag[i] = j; 1718 break; 1719 } 1720 } 1721 } 1722 PetscFunctionReturn(0); 1723 } 1724 1725 PetscErrorCode MatShift_SeqAIJ(Mat A,PetscScalar v) 1726 { 1727 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1728 const PetscInt *diag = (const PetscInt*)a->diag; 1729 const PetscInt *ii = (const PetscInt*) a->i; 1730 PetscInt i,*mdiag = NULL; 1731 PetscInt cnt = 0; /* how many diagonals are missing */ 1732 1733 PetscFunctionBegin; 1734 if (!A->preallocated || !a->nz) { 1735 PetscCall(MatSeqAIJSetPreallocation(A,1,NULL)); 1736 PetscCall(MatShift_Basic(A,v)); 1737 PetscFunctionReturn(0); 1738 } 1739 1740 if (a->diagonaldense) { 1741 cnt = 0; 1742 } else { 1743 PetscCall(PetscCalloc1(A->rmap->n,&mdiag)); 1744 for (i=0; i<A->rmap->n; i++) { 1745 if (i < A->cmap->n && diag[i] >= ii[i+1]) { /* 'out of range' rows never have diagonals */ 1746 cnt++; 1747 mdiag[i] = 1; 1748 } 1749 } 1750 } 1751 if (!cnt) { 1752 PetscCall(MatShift_Basic(A,v)); 1753 } else { 1754 PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */ 1755 PetscInt *oldj = a->j, *oldi = a->i; 1756 PetscBool singlemalloc = a->singlemalloc,free_a = a->free_a,free_ij = a->free_ij; 1757 1758 a->a = NULL; 1759 a->j = NULL; 1760 a->i = NULL; 1761 /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */ 1762 for (i=0; i<PetscMin(A->rmap->n,A->cmap->n); i++) { 1763 a->imax[i] += mdiag[i]; 1764 } 1765 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax)); 1766 1767 /* copy old values into new matrix data structure */ 1768 for (i=0; i<A->rmap->n; i++) { 1769 PetscCall(MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES)); 1770 if (i < A->cmap->n) { 1771 PetscCall(MatSetValue(A,i,i,v,ADD_VALUES)); 1772 } 1773 } 1774 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1775 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1776 if (singlemalloc) { 1777 PetscCall(PetscFree3(olda,oldj,oldi)); 1778 } else { 1779 if (free_a) PetscCall(PetscFree(olda)); 1780 if (free_ij) PetscCall(PetscFree(oldj)); 1781 if (free_ij) PetscCall(PetscFree(oldi)); 1782 } 1783 } 1784 PetscCall(PetscFree(mdiag)); 1785 a->diagonaldense = PETSC_TRUE; 1786 PetscFunctionReturn(0); 1787 } 1788 1789 /* 1790 Checks for missing diagonals 1791 */ 1792 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1793 { 1794 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1795 PetscInt *diag,*ii = a->i,i; 1796 1797 PetscFunctionBegin; 1798 *missing = PETSC_FALSE; 1799 if (A->rmap->n > 0 && !ii) { 1800 *missing = PETSC_TRUE; 1801 if (d) *d = 0; 1802 PetscCall(PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n")); 1803 } else { 1804 PetscInt n; 1805 n = PetscMin(A->rmap->n, A->cmap->n); 1806 diag = a->diag; 1807 for (i=0; i<n; i++) { 1808 if (diag[i] >= ii[i+1]) { 1809 *missing = PETSC_TRUE; 1810 if (d) *d = i; 1811 PetscCall(PetscInfo(A,"Matrix is missing diagonal number %" PetscInt_FMT "\n",i)); 1812 break; 1813 } 1814 } 1815 } 1816 PetscFunctionReturn(0); 1817 } 1818 1819 #include <petscblaslapack.h> 1820 #include <petsc/private/kernels/blockinvert.h> 1821 1822 /* 1823 Note that values is allocated externally by the PC and then passed into this routine 1824 */ 1825 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag) 1826 { 1827 PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots; 1828 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1829 const PetscReal shift = 0.0; 1830 PetscInt ipvt[5]; 1831 PetscScalar work[25],*v_work; 1832 1833 PetscFunctionBegin; 1834 allowzeropivot = PetscNot(A->erroriffailure); 1835 for (i=0; i<nblocks; i++) ncnt += bsizes[i]; 1836 PetscCheck(ncnt == n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Total blocksizes %" PetscInt_FMT " doesn't match number matrix rows %" PetscInt_FMT,ncnt,n); 1837 for (i=0; i<nblocks; i++) { 1838 bsizemax = PetscMax(bsizemax,bsizes[i]); 1839 } 1840 PetscCall(PetscMalloc1(bsizemax,&indx)); 1841 if (bsizemax > 7) { 1842 PetscCall(PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots)); 1843 } 1844 ncnt = 0; 1845 for (i=0; i<nblocks; i++) { 1846 for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j; 1847 PetscCall(MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag)); 1848 switch (bsizes[i]) { 1849 case 1: 1850 *diag = 1.0/(*diag); 1851 break; 1852 case 2: 1853 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 1854 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1855 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 1856 break; 1857 case 3: 1858 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 1859 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1860 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 1861 break; 1862 case 4: 1863 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 1864 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1865 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 1866 break; 1867 case 5: 1868 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 1869 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1870 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 1871 break; 1872 case 6: 1873 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 1874 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1875 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 1876 break; 1877 case 7: 1878 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 1879 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1880 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 1881 break; 1882 default: 1883 PetscCall(PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 1884 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1885 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bsizes[i])); 1886 } 1887 ncnt += bsizes[i]; 1888 diag += bsizes[i]*bsizes[i]; 1889 } 1890 if (bsizemax > 7) { 1891 PetscCall(PetscFree2(v_work,v_pivots)); 1892 } 1893 PetscCall(PetscFree(indx)); 1894 PetscFunctionReturn(0); 1895 } 1896 1897 /* 1898 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1899 */ 1900 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1901 { 1902 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1903 PetscInt i,*diag,m = A->rmap->n; 1904 const MatScalar *v; 1905 PetscScalar *idiag,*mdiag; 1906 1907 PetscFunctionBegin; 1908 if (a->idiagvalid) PetscFunctionReturn(0); 1909 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1910 diag = a->diag; 1911 if (!a->idiag) { 1912 PetscCall(PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work)); 1913 PetscCall(PetscLogObjectMemory((PetscObject)A,3*m*sizeof(PetscScalar))); 1914 } 1915 1916 mdiag = a->mdiag; 1917 idiag = a->idiag; 1918 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 1919 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1920 for (i=0; i<m; i++) { 1921 mdiag[i] = v[diag[i]]; 1922 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1923 if (PetscRealPart(fshift)) { 1924 PetscCall(PetscInfo(A,"Zero diagonal on row %" PetscInt_FMT "\n",i)); 1925 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1926 A->factorerror_zeropivot_value = 0.0; 1927 A->factorerror_zeropivot_row = i; 1928 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %" PetscInt_FMT,i); 1929 } 1930 idiag[i] = 1.0/v[diag[i]]; 1931 } 1932 PetscCall(PetscLogFlops(m)); 1933 } else { 1934 for (i=0; i<m; i++) { 1935 mdiag[i] = v[diag[i]]; 1936 idiag[i] = omega/(fshift + v[diag[i]]); 1937 } 1938 PetscCall(PetscLogFlops(2.0*m)); 1939 } 1940 a->idiagvalid = PETSC_TRUE; 1941 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 1942 PetscFunctionReturn(0); 1943 } 1944 1945 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1946 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1947 { 1948 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1949 PetscScalar *x,d,sum,*t,scale; 1950 const MatScalar *v,*idiag=NULL,*mdiag,*aa; 1951 const PetscScalar *b, *bs,*xb, *ts; 1952 PetscInt n,m = A->rmap->n,i; 1953 const PetscInt *idx,*diag; 1954 1955 PetscFunctionBegin; 1956 if (a->inode.use && a->inode.checked && omega == 1.0 && fshift == 0.0) { 1957 PetscCall(MatSOR_SeqAIJ_Inode(A,bb,omega,flag,fshift,its,lits,xx)); 1958 PetscFunctionReturn(0); 1959 } 1960 its = its*lits; 1961 1962 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1963 if (!a->idiagvalid) PetscCall(MatInvertDiagonal_SeqAIJ(A,omega,fshift)); 1964 a->fshift = fshift; 1965 a->omega = omega; 1966 1967 diag = a->diag; 1968 t = a->ssor_work; 1969 idiag = a->idiag; 1970 mdiag = a->mdiag; 1971 1972 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1973 PetscCall(VecGetArray(xx,&x)); 1974 PetscCall(VecGetArrayRead(bb,&b)); 1975 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1976 if (flag == SOR_APPLY_UPPER) { 1977 /* apply (U + D/omega) to the vector */ 1978 bs = b; 1979 for (i=0; i<m; i++) { 1980 d = fshift + mdiag[i]; 1981 n = a->i[i+1] - diag[i] - 1; 1982 idx = a->j + diag[i] + 1; 1983 v = aa + diag[i] + 1; 1984 sum = b[i]*d/omega; 1985 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1986 x[i] = sum; 1987 } 1988 PetscCall(VecRestoreArray(xx,&x)); 1989 PetscCall(VecRestoreArrayRead(bb,&b)); 1990 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1991 PetscCall(PetscLogFlops(a->nz)); 1992 PetscFunctionReturn(0); 1993 } 1994 1995 PetscCheck(flag != SOR_APPLY_LOWER,PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1996 if (flag & SOR_EISENSTAT) { 1997 /* Let A = L + U + D; where L is lower triangular, 1998 U is upper triangular, E = D/omega; This routine applies 1999 2000 (L + E)^{-1} A (U + E)^{-1} 2001 2002 to a vector efficiently using Eisenstat's trick. 2003 */ 2004 scale = (2.0/omega) - 1.0; 2005 2006 /* x = (E + U)^{-1} b */ 2007 for (i=m-1; i>=0; i--) { 2008 n = a->i[i+1] - diag[i] - 1; 2009 idx = a->j + diag[i] + 1; 2010 v = aa + diag[i] + 1; 2011 sum = b[i]; 2012 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2013 x[i] = sum*idiag[i]; 2014 } 2015 2016 /* t = b - (2*E - D)x */ 2017 v = aa; 2018 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 2019 2020 /* t = (E + L)^{-1}t */ 2021 ts = t; 2022 diag = a->diag; 2023 for (i=0; i<m; i++) { 2024 n = diag[i] - a->i[i]; 2025 idx = a->j + a->i[i]; 2026 v = aa + a->i[i]; 2027 sum = t[i]; 2028 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 2029 t[i] = sum*idiag[i]; 2030 /* x = x + t */ 2031 x[i] += t[i]; 2032 } 2033 2034 PetscCall(PetscLogFlops(6.0*m-1 + 2.0*a->nz)); 2035 PetscCall(VecRestoreArray(xx,&x)); 2036 PetscCall(VecRestoreArrayRead(bb,&b)); 2037 PetscFunctionReturn(0); 2038 } 2039 if (flag & SOR_ZERO_INITIAL_GUESS) { 2040 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2041 for (i=0; i<m; i++) { 2042 n = diag[i] - a->i[i]; 2043 idx = a->j + a->i[i]; 2044 v = aa + a->i[i]; 2045 sum = b[i]; 2046 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2047 t[i] = sum; 2048 x[i] = sum*idiag[i]; 2049 } 2050 xb = t; 2051 PetscCall(PetscLogFlops(a->nz)); 2052 } else xb = b; 2053 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2054 for (i=m-1; i>=0; i--) { 2055 n = a->i[i+1] - diag[i] - 1; 2056 idx = a->j + diag[i] + 1; 2057 v = aa + diag[i] + 1; 2058 sum = xb[i]; 2059 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2060 if (xb == b) { 2061 x[i] = sum*idiag[i]; 2062 } else { 2063 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2064 } 2065 } 2066 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2067 } 2068 its--; 2069 } 2070 while (its--) { 2071 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2072 for (i=0; i<m; i++) { 2073 /* lower */ 2074 n = diag[i] - a->i[i]; 2075 idx = a->j + a->i[i]; 2076 v = aa + a->i[i]; 2077 sum = b[i]; 2078 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2079 t[i] = sum; /* save application of the lower-triangular part */ 2080 /* upper */ 2081 n = a->i[i+1] - diag[i] - 1; 2082 idx = a->j + diag[i] + 1; 2083 v = aa + diag[i] + 1; 2084 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2085 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2086 } 2087 xb = t; 2088 PetscCall(PetscLogFlops(2.0*a->nz)); 2089 } else xb = b; 2090 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2091 for (i=m-1; i>=0; i--) { 2092 sum = xb[i]; 2093 if (xb == b) { 2094 /* whole matrix (no checkpointing available) */ 2095 n = a->i[i+1] - a->i[i]; 2096 idx = a->j + a->i[i]; 2097 v = aa + a->i[i]; 2098 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2099 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 2100 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 2101 n = a->i[i+1] - diag[i] - 1; 2102 idx = a->j + diag[i] + 1; 2103 v = aa + diag[i] + 1; 2104 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2105 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2106 } 2107 } 2108 if (xb == b) { 2109 PetscCall(PetscLogFlops(2.0*a->nz)); 2110 } else { 2111 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2112 } 2113 } 2114 } 2115 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2116 PetscCall(VecRestoreArray(xx,&x)); 2117 PetscCall(VecRestoreArrayRead(bb,&b)); 2118 PetscFunctionReturn(0); 2119 } 2120 2121 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 2122 { 2123 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2124 2125 PetscFunctionBegin; 2126 info->block_size = 1.0; 2127 info->nz_allocated = a->maxnz; 2128 info->nz_used = a->nz; 2129 info->nz_unneeded = (a->maxnz - a->nz); 2130 info->assemblies = A->num_ass; 2131 info->mallocs = A->info.mallocs; 2132 info->memory = ((PetscObject)A)->mem; 2133 if (A->factortype) { 2134 info->fill_ratio_given = A->info.fill_ratio_given; 2135 info->fill_ratio_needed = A->info.fill_ratio_needed; 2136 info->factor_mallocs = A->info.factor_mallocs; 2137 } else { 2138 info->fill_ratio_given = 0; 2139 info->fill_ratio_needed = 0; 2140 info->factor_mallocs = 0; 2141 } 2142 PetscFunctionReturn(0); 2143 } 2144 2145 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2146 { 2147 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2148 PetscInt i,m = A->rmap->n - 1; 2149 const PetscScalar *xx; 2150 PetscScalar *bb,*aa; 2151 PetscInt d = 0; 2152 2153 PetscFunctionBegin; 2154 if (x && b) { 2155 PetscCall(VecGetArrayRead(x,&xx)); 2156 PetscCall(VecGetArray(b,&bb)); 2157 for (i=0; i<N; i++) { 2158 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2159 if (rows[i] >= A->cmap->n) continue; 2160 bb[rows[i]] = diag*xx[rows[i]]; 2161 } 2162 PetscCall(VecRestoreArrayRead(x,&xx)); 2163 PetscCall(VecRestoreArray(b,&bb)); 2164 } 2165 2166 PetscCall(MatSeqAIJGetArray(A,&aa)); 2167 if (a->keepnonzeropattern) { 2168 for (i=0; i<N; i++) { 2169 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2170 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2171 } 2172 if (diag != 0.0) { 2173 for (i=0; i<N; i++) { 2174 d = rows[i]; 2175 if (rows[i] >= A->cmap->n) continue; 2176 PetscCheck(a->diag[d] < a->i[d+1],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in the zeroed row %" PetscInt_FMT,d); 2177 } 2178 for (i=0; i<N; i++) { 2179 if (rows[i] >= A->cmap->n) continue; 2180 aa[a->diag[rows[i]]] = diag; 2181 } 2182 } 2183 } else { 2184 if (diag != 0.0) { 2185 for (i=0; i<N; i++) { 2186 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2187 if (a->ilen[rows[i]] > 0) { 2188 if (rows[i] >= A->cmap->n) { 2189 a->ilen[rows[i]] = 0; 2190 } else { 2191 a->ilen[rows[i]] = 1; 2192 aa[a->i[rows[i]]] = diag; 2193 a->j[a->i[rows[i]]] = rows[i]; 2194 } 2195 } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */ 2196 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2197 } 2198 } 2199 } else { 2200 for (i=0; i<N; i++) { 2201 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2202 a->ilen[rows[i]] = 0; 2203 } 2204 } 2205 A->nonzerostate++; 2206 } 2207 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2208 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2209 PetscFunctionReturn(0); 2210 } 2211 2212 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2213 { 2214 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2215 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2216 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2217 const PetscScalar *xx; 2218 PetscScalar *bb,*aa; 2219 2220 PetscFunctionBegin; 2221 if (!N) PetscFunctionReturn(0); 2222 PetscCall(MatSeqAIJGetArray(A,&aa)); 2223 if (x && b) { 2224 PetscCall(VecGetArrayRead(x,&xx)); 2225 PetscCall(VecGetArray(b,&bb)); 2226 vecs = PETSC_TRUE; 2227 } 2228 PetscCall(PetscCalloc1(A->rmap->n,&zeroed)); 2229 for (i=0; i<N; i++) { 2230 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2231 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2232 2233 zeroed[rows[i]] = PETSC_TRUE; 2234 } 2235 for (i=0; i<A->rmap->n; i++) { 2236 if (!zeroed[i]) { 2237 for (j=a->i[i]; j<a->i[i+1]; j++) { 2238 if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) { 2239 if (vecs) bb[i] -= aa[j]*xx[a->j[j]]; 2240 aa[j] = 0.0; 2241 } 2242 } 2243 } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i]; 2244 } 2245 if (x && b) { 2246 PetscCall(VecRestoreArrayRead(x,&xx)); 2247 PetscCall(VecRestoreArray(b,&bb)); 2248 } 2249 PetscCall(PetscFree(zeroed)); 2250 if (diag != 0.0) { 2251 PetscCall(MatMissingDiagonal_SeqAIJ(A,&missing,&d)); 2252 if (missing) { 2253 for (i=0; i<N; i++) { 2254 if (rows[i] >= A->cmap->N) continue; 2255 PetscCheck(!a->nonew || rows[i] < d,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %" PetscInt_FMT " (%" PetscInt_FMT ")",d,rows[i]); 2256 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2257 } 2258 } else { 2259 for (i=0; i<N; i++) { 2260 aa[a->diag[rows[i]]] = diag; 2261 } 2262 } 2263 } 2264 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2265 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2266 PetscFunctionReturn(0); 2267 } 2268 2269 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2270 { 2271 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2272 const PetscScalar *aa; 2273 PetscInt *itmp; 2274 2275 PetscFunctionBegin; 2276 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2277 *nz = a->i[row+1] - a->i[row]; 2278 if (v) *v = (PetscScalar*)(aa + a->i[row]); 2279 if (idx) { 2280 itmp = a->j + a->i[row]; 2281 if (*nz) *idx = itmp; 2282 else *idx = NULL; 2283 } 2284 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2285 PetscFunctionReturn(0); 2286 } 2287 2288 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2289 { 2290 PetscFunctionBegin; 2291 if (nz) *nz = 0; 2292 if (idx) *idx = NULL; 2293 if (v) *v = NULL; 2294 PetscFunctionReturn(0); 2295 } 2296 2297 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2298 { 2299 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2300 const MatScalar *v; 2301 PetscReal sum = 0.0; 2302 PetscInt i,j; 2303 2304 PetscFunctionBegin; 2305 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 2306 if (type == NORM_FROBENIUS) { 2307 #if defined(PETSC_USE_REAL___FP16) 2308 PetscBLASInt one = 1,nz = a->nz; 2309 PetscCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&nz,v,&one)); 2310 #else 2311 for (i=0; i<a->nz; i++) { 2312 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2313 } 2314 *nrm = PetscSqrtReal(sum); 2315 #endif 2316 PetscCall(PetscLogFlops(2.0*a->nz)); 2317 } else if (type == NORM_1) { 2318 PetscReal *tmp; 2319 PetscInt *jj = a->j; 2320 PetscCall(PetscCalloc1(A->cmap->n+1,&tmp)); 2321 *nrm = 0.0; 2322 for (j=0; j<a->nz; j++) { 2323 tmp[*jj++] += PetscAbsScalar(*v); v++; 2324 } 2325 for (j=0; j<A->cmap->n; j++) { 2326 if (tmp[j] > *nrm) *nrm = tmp[j]; 2327 } 2328 PetscCall(PetscFree(tmp)); 2329 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2330 } else if (type == NORM_INFINITY) { 2331 *nrm = 0.0; 2332 for (j=0; j<A->rmap->n; j++) { 2333 const PetscScalar *v2 = v + a->i[j]; 2334 sum = 0.0; 2335 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2336 sum += PetscAbsScalar(*v2); v2++; 2337 } 2338 if (sum > *nrm) *nrm = sum; 2339 } 2340 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2341 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2342 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 2343 PetscFunctionReturn(0); 2344 } 2345 2346 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2347 { 2348 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2349 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2350 const MatScalar *va,*vb; 2351 PetscInt ma,na,mb,nb, i; 2352 2353 PetscFunctionBegin; 2354 PetscCall(MatGetSize(A,&ma,&na)); 2355 PetscCall(MatGetSize(B,&mb,&nb)); 2356 if (ma!=nb || na!=mb) { 2357 *f = PETSC_FALSE; 2358 PetscFunctionReturn(0); 2359 } 2360 PetscCall(MatSeqAIJGetArrayRead(A,&va)); 2361 PetscCall(MatSeqAIJGetArrayRead(B,&vb)); 2362 aii = aij->i; bii = bij->i; 2363 adx = aij->j; bdx = bij->j; 2364 PetscCall(PetscMalloc1(ma,&aptr)); 2365 PetscCall(PetscMalloc1(mb,&bptr)); 2366 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2367 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2368 2369 *f = PETSC_TRUE; 2370 for (i=0; i<ma; i++) { 2371 while (aptr[i]<aii[i+1]) { 2372 PetscInt idc,idr; 2373 PetscScalar vc,vr; 2374 /* column/row index/value */ 2375 idc = adx[aptr[i]]; 2376 idr = bdx[bptr[idc]]; 2377 vc = va[aptr[i]]; 2378 vr = vb[bptr[idc]]; 2379 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2380 *f = PETSC_FALSE; 2381 goto done; 2382 } else { 2383 aptr[i]++; 2384 if (B || i!=idc) bptr[idc]++; 2385 } 2386 } 2387 } 2388 done: 2389 PetscCall(PetscFree(aptr)); 2390 PetscCall(PetscFree(bptr)); 2391 PetscCall(MatSeqAIJRestoreArrayRead(A,&va)); 2392 PetscCall(MatSeqAIJRestoreArrayRead(B,&vb)); 2393 PetscFunctionReturn(0); 2394 } 2395 2396 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2397 { 2398 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2399 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2400 MatScalar *va,*vb; 2401 PetscInt ma,na,mb,nb, i; 2402 2403 PetscFunctionBegin; 2404 PetscCall(MatGetSize(A,&ma,&na)); 2405 PetscCall(MatGetSize(B,&mb,&nb)); 2406 if (ma!=nb || na!=mb) { 2407 *f = PETSC_FALSE; 2408 PetscFunctionReturn(0); 2409 } 2410 aii = aij->i; bii = bij->i; 2411 adx = aij->j; bdx = bij->j; 2412 va = aij->a; vb = bij->a; 2413 PetscCall(PetscMalloc1(ma,&aptr)); 2414 PetscCall(PetscMalloc1(mb,&bptr)); 2415 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2416 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2417 2418 *f = PETSC_TRUE; 2419 for (i=0; i<ma; i++) { 2420 while (aptr[i]<aii[i+1]) { 2421 PetscInt idc,idr; 2422 PetscScalar vc,vr; 2423 /* column/row index/value */ 2424 idc = adx[aptr[i]]; 2425 idr = bdx[bptr[idc]]; 2426 vc = va[aptr[i]]; 2427 vr = vb[bptr[idc]]; 2428 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2429 *f = PETSC_FALSE; 2430 goto done; 2431 } else { 2432 aptr[i]++; 2433 if (B || i!=idc) bptr[idc]++; 2434 } 2435 } 2436 } 2437 done: 2438 PetscCall(PetscFree(aptr)); 2439 PetscCall(PetscFree(bptr)); 2440 PetscFunctionReturn(0); 2441 } 2442 2443 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2444 { 2445 PetscFunctionBegin; 2446 PetscCall(MatIsTranspose_SeqAIJ(A,A,tol,f)); 2447 PetscFunctionReturn(0); 2448 } 2449 2450 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2451 { 2452 PetscFunctionBegin; 2453 PetscCall(MatIsHermitianTranspose_SeqAIJ(A,A,tol,f)); 2454 PetscFunctionReturn(0); 2455 } 2456 2457 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2458 { 2459 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2460 const PetscScalar *l,*r; 2461 PetscScalar x; 2462 MatScalar *v; 2463 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2464 const PetscInt *jj; 2465 2466 PetscFunctionBegin; 2467 if (ll) { 2468 /* The local size is used so that VecMPI can be passed to this routine 2469 by MatDiagonalScale_MPIAIJ */ 2470 PetscCall(VecGetLocalSize(ll,&m)); 2471 PetscCheck(m == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2472 PetscCall(VecGetArrayRead(ll,&l)); 2473 PetscCall(MatSeqAIJGetArray(A,&v)); 2474 for (i=0; i<m; i++) { 2475 x = l[i]; 2476 M = a->i[i+1] - a->i[i]; 2477 for (j=0; j<M; j++) (*v++) *= x; 2478 } 2479 PetscCall(VecRestoreArrayRead(ll,&l)); 2480 PetscCall(PetscLogFlops(nz)); 2481 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2482 } 2483 if (rr) { 2484 PetscCall(VecGetLocalSize(rr,&n)); 2485 PetscCheck(n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2486 PetscCall(VecGetArrayRead(rr,&r)); 2487 PetscCall(MatSeqAIJGetArray(A,&v)); 2488 jj = a->j; 2489 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2490 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2491 PetscCall(VecRestoreArrayRead(rr,&r)); 2492 PetscCall(PetscLogFlops(nz)); 2493 } 2494 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 2495 PetscFunctionReturn(0); 2496 } 2497 2498 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2499 { 2500 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2501 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2502 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2503 const PetscInt *irow,*icol; 2504 const PetscScalar *aa; 2505 PetscInt nrows,ncols; 2506 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2507 MatScalar *a_new,*mat_a; 2508 Mat C; 2509 PetscBool stride; 2510 2511 PetscFunctionBegin; 2512 PetscCall(ISGetIndices(isrow,&irow)); 2513 PetscCall(ISGetLocalSize(isrow,&nrows)); 2514 PetscCall(ISGetLocalSize(iscol,&ncols)); 2515 2516 PetscCall(PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride)); 2517 if (stride) { 2518 PetscCall(ISStrideGetInfo(iscol,&first,&step)); 2519 } else { 2520 first = 0; 2521 step = 0; 2522 } 2523 if (stride && step == 1) { 2524 /* special case of contiguous rows */ 2525 PetscCall(PetscMalloc2(nrows,&lens,nrows,&starts)); 2526 /* loop over new rows determining lens and starting points */ 2527 for (i=0; i<nrows; i++) { 2528 kstart = ai[irow[i]]; 2529 kend = kstart + ailen[irow[i]]; 2530 starts[i] = kstart; 2531 for (k=kstart; k<kend; k++) { 2532 if (aj[k] >= first) { 2533 starts[i] = k; 2534 break; 2535 } 2536 } 2537 sum = 0; 2538 while (k < kend) { 2539 if (aj[k++] >= first+ncols) break; 2540 sum++; 2541 } 2542 lens[i] = sum; 2543 } 2544 /* create submatrix */ 2545 if (scall == MAT_REUSE_MATRIX) { 2546 PetscInt n_cols,n_rows; 2547 PetscCall(MatGetSize(*B,&n_rows,&n_cols)); 2548 PetscCheck(n_rows == nrows && n_cols == ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2549 PetscCall(MatZeroEntries(*B)); 2550 C = *B; 2551 } else { 2552 PetscInt rbs,cbs; 2553 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2554 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2555 PetscCall(ISGetBlockSize(isrow,&rbs)); 2556 PetscCall(ISGetBlockSize(iscol,&cbs)); 2557 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2558 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2559 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2560 } 2561 c = (Mat_SeqAIJ*)C->data; 2562 2563 /* loop over rows inserting into submatrix */ 2564 a_new = c->a; 2565 j_new = c->j; 2566 i_new = c->i; 2567 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2568 for (i=0; i<nrows; i++) { 2569 ii = starts[i]; 2570 lensi = lens[i]; 2571 for (k=0; k<lensi; k++) { 2572 *j_new++ = aj[ii+k] - first; 2573 } 2574 PetscCall(PetscArraycpy(a_new,aa + starts[i],lensi)); 2575 a_new += lensi; 2576 i_new[i+1] = i_new[i] + lensi; 2577 c->ilen[i] = lensi; 2578 } 2579 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2580 PetscCall(PetscFree2(lens,starts)); 2581 } else { 2582 PetscCall(ISGetIndices(iscol,&icol)); 2583 PetscCall(PetscCalloc1(oldcols,&smap)); 2584 PetscCall(PetscMalloc1(1+nrows,&lens)); 2585 for (i=0; i<ncols; i++) { 2586 PetscCheck(icol[i] < oldcols,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%" PetscInt_FMT "] %" PetscInt_FMT " >= A->cmap->n %" PetscInt_FMT,i,icol[i],oldcols); 2587 smap[icol[i]] = i+1; 2588 } 2589 2590 /* determine lens of each row */ 2591 for (i=0; i<nrows; i++) { 2592 kstart = ai[irow[i]]; 2593 kend = kstart + a->ilen[irow[i]]; 2594 lens[i] = 0; 2595 for (k=kstart; k<kend; k++) { 2596 if (smap[aj[k]]) { 2597 lens[i]++; 2598 } 2599 } 2600 } 2601 /* Create and fill new matrix */ 2602 if (scall == MAT_REUSE_MATRIX) { 2603 PetscBool equal; 2604 2605 c = (Mat_SeqAIJ*)((*B)->data); 2606 PetscCheck((*B)->rmap->n == nrows && (*B)->cmap->n == ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2607 PetscCall(PetscArraycmp(c->ilen,lens,(*B)->rmap->n,&equal)); 2608 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2609 PetscCall(PetscArrayzero(c->ilen,(*B)->rmap->n)); 2610 C = *B; 2611 } else { 2612 PetscInt rbs,cbs; 2613 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2614 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2615 PetscCall(ISGetBlockSize(isrow,&rbs)); 2616 PetscCall(ISGetBlockSize(iscol,&cbs)); 2617 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2618 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2619 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2620 } 2621 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2622 c = (Mat_SeqAIJ*)(C->data); 2623 for (i=0; i<nrows; i++) { 2624 row = irow[i]; 2625 kstart = ai[row]; 2626 kend = kstart + a->ilen[row]; 2627 mat_i = c->i[i]; 2628 mat_j = c->j + mat_i; 2629 mat_a = c->a + mat_i; 2630 mat_ilen = c->ilen + i; 2631 for (k=kstart; k<kend; k++) { 2632 if ((tcol=smap[a->j[k]])) { 2633 *mat_j++ = tcol - 1; 2634 *mat_a++ = aa[k]; 2635 (*mat_ilen)++; 2636 2637 } 2638 } 2639 } 2640 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2641 /* Free work space */ 2642 PetscCall(ISRestoreIndices(iscol,&icol)); 2643 PetscCall(PetscFree(smap)); 2644 PetscCall(PetscFree(lens)); 2645 /* sort */ 2646 for (i = 0; i < nrows; i++) { 2647 PetscInt ilen; 2648 2649 mat_i = c->i[i]; 2650 mat_j = c->j + mat_i; 2651 mat_a = c->a + mat_i; 2652 ilen = c->ilen[i]; 2653 PetscCall(PetscSortIntWithScalarArray(ilen,mat_j,mat_a)); 2654 } 2655 } 2656 #if defined(PETSC_HAVE_DEVICE) 2657 PetscCall(MatBindToCPU(C,A->boundtocpu)); 2658 #endif 2659 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 2660 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 2661 2662 PetscCall(ISRestoreIndices(isrow,&irow)); 2663 *B = C; 2664 PetscFunctionReturn(0); 2665 } 2666 2667 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2668 { 2669 Mat B; 2670 2671 PetscFunctionBegin; 2672 if (scall == MAT_INITIAL_MATRIX) { 2673 PetscCall(MatCreate(subComm,&B)); 2674 PetscCall(MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n)); 2675 PetscCall(MatSetBlockSizesFromMats(B,mat,mat)); 2676 PetscCall(MatSetType(B,MATSEQAIJ)); 2677 PetscCall(MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE)); 2678 *subMat = B; 2679 } else { 2680 PetscCall(MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN)); 2681 } 2682 PetscFunctionReturn(0); 2683 } 2684 2685 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2686 { 2687 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2688 Mat outA; 2689 PetscBool row_identity,col_identity; 2690 2691 PetscFunctionBegin; 2692 PetscCheck(info->levels == 0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2693 2694 PetscCall(ISIdentity(row,&row_identity)); 2695 PetscCall(ISIdentity(col,&col_identity)); 2696 2697 outA = inA; 2698 outA->factortype = MAT_FACTOR_LU; 2699 PetscCall(PetscFree(inA->solvertype)); 2700 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype)); 2701 2702 PetscCall(PetscObjectReference((PetscObject)row)); 2703 PetscCall(ISDestroy(&a->row)); 2704 2705 a->row = row; 2706 2707 PetscCall(PetscObjectReference((PetscObject)col)); 2708 PetscCall(ISDestroy(&a->col)); 2709 2710 a->col = col; 2711 2712 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2713 PetscCall(ISDestroy(&a->icol)); 2714 PetscCall(ISInvertPermutation(col,PETSC_DECIDE,&a->icol)); 2715 PetscCall(PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol)); 2716 2717 if (!a->solve_work) { /* this matrix may have been factored before */ 2718 PetscCall(PetscMalloc1(inA->rmap->n+1,&a->solve_work)); 2719 PetscCall(PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar))); 2720 } 2721 2722 PetscCall(MatMarkDiagonal_SeqAIJ(inA)); 2723 if (row_identity && col_identity) { 2724 PetscCall(MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info)); 2725 } else { 2726 PetscCall(MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info)); 2727 } 2728 PetscFunctionReturn(0); 2729 } 2730 2731 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2732 { 2733 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2734 PetscScalar *v; 2735 PetscBLASInt one = 1,bnz; 2736 2737 PetscFunctionBegin; 2738 PetscCall(MatSeqAIJGetArray(inA,&v)); 2739 PetscCall(PetscBLASIntCast(a->nz,&bnz)); 2740 PetscCallBLAS("BLASscal",BLASscal_(&bnz,&alpha,v,&one)); 2741 PetscCall(PetscLogFlops(a->nz)); 2742 PetscCall(MatSeqAIJRestoreArray(inA,&v)); 2743 PetscCall(MatSeqAIJInvalidateDiagonal(inA)); 2744 PetscFunctionReturn(0); 2745 } 2746 2747 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2748 { 2749 PetscInt i; 2750 2751 PetscFunctionBegin; 2752 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2753 PetscCall(PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr)); 2754 2755 for (i=0; i<submatj->nrqr; ++i) { 2756 PetscCall(PetscFree(submatj->sbuf2[i])); 2757 } 2758 PetscCall(PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1)); 2759 2760 if (submatj->rbuf1) { 2761 PetscCall(PetscFree(submatj->rbuf1[0])); 2762 PetscCall(PetscFree(submatj->rbuf1)); 2763 } 2764 2765 for (i=0; i<submatj->nrqs; ++i) { 2766 PetscCall(PetscFree(submatj->rbuf3[i])); 2767 } 2768 PetscCall(PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3)); 2769 PetscCall(PetscFree(submatj->pa)); 2770 } 2771 2772 #if defined(PETSC_USE_CTABLE) 2773 PetscCall(PetscTableDestroy((PetscTable*)&submatj->rmap)); 2774 if (submatj->cmap_loc) PetscCall(PetscFree(submatj->cmap_loc)); 2775 PetscCall(PetscFree(submatj->rmap_loc)); 2776 #else 2777 PetscCall(PetscFree(submatj->rmap)); 2778 #endif 2779 2780 if (!submatj->allcolumns) { 2781 #if defined(PETSC_USE_CTABLE) 2782 PetscCall(PetscTableDestroy((PetscTable*)&submatj->cmap)); 2783 #else 2784 PetscCall(PetscFree(submatj->cmap)); 2785 #endif 2786 } 2787 PetscCall(PetscFree(submatj->row2proc)); 2788 2789 PetscCall(PetscFree(submatj)); 2790 PetscFunctionReturn(0); 2791 } 2792 2793 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2794 { 2795 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2796 Mat_SubSppt *submatj = c->submatis1; 2797 2798 PetscFunctionBegin; 2799 PetscCall((*submatj->destroy)(C)); 2800 PetscCall(MatDestroySubMatrix_Private(submatj)); 2801 PetscFunctionReturn(0); 2802 } 2803 2804 /* Note this has code duplication with MatDestroySubMatrices_SeqBAIJ() */ 2805 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2806 { 2807 PetscInt i; 2808 Mat C; 2809 Mat_SeqAIJ *c; 2810 Mat_SubSppt *submatj; 2811 2812 PetscFunctionBegin; 2813 for (i=0; i<n; i++) { 2814 C = (*mat)[i]; 2815 c = (Mat_SeqAIJ*)C->data; 2816 submatj = c->submatis1; 2817 if (submatj) { 2818 if (--((PetscObject)C)->refct <= 0) { 2819 PetscCall(PetscFree(C->factorprefix)); 2820 PetscCall((*submatj->destroy)(C)); 2821 PetscCall(MatDestroySubMatrix_Private(submatj)); 2822 PetscCall(PetscFree(C->defaultvectype)); 2823 PetscCall(PetscLayoutDestroy(&C->rmap)); 2824 PetscCall(PetscLayoutDestroy(&C->cmap)); 2825 PetscCall(PetscHeaderDestroy(&C)); 2826 } 2827 } else { 2828 PetscCall(MatDestroy(&C)); 2829 } 2830 } 2831 2832 /* Destroy Dummy submatrices created for reuse */ 2833 PetscCall(MatDestroySubMatrices_Dummy(n,mat)); 2834 2835 PetscCall(PetscFree(*mat)); 2836 PetscFunctionReturn(0); 2837 } 2838 2839 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2840 { 2841 PetscInt i; 2842 2843 PetscFunctionBegin; 2844 if (scall == MAT_INITIAL_MATRIX) { 2845 PetscCall(PetscCalloc1(n+1,B)); 2846 } 2847 2848 for (i=0; i<n; i++) { 2849 PetscCall(MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i])); 2850 } 2851 PetscFunctionReturn(0); 2852 } 2853 2854 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2855 { 2856 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2857 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2858 const PetscInt *idx; 2859 PetscInt start,end,*ai,*aj; 2860 PetscBT table; 2861 2862 PetscFunctionBegin; 2863 m = A->rmap->n; 2864 ai = a->i; 2865 aj = a->j; 2866 2867 PetscCheck(ov >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2868 2869 PetscCall(PetscMalloc1(m+1,&nidx)); 2870 PetscCall(PetscBTCreate(m,&table)); 2871 2872 for (i=0; i<is_max; i++) { 2873 /* Initialize the two local arrays */ 2874 isz = 0; 2875 PetscCall(PetscBTMemzero(m,table)); 2876 2877 /* Extract the indices, assume there can be duplicate entries */ 2878 PetscCall(ISGetIndices(is[i],&idx)); 2879 PetscCall(ISGetLocalSize(is[i],&n)); 2880 2881 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2882 for (j=0; j<n; ++j) { 2883 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2884 } 2885 PetscCall(ISRestoreIndices(is[i],&idx)); 2886 PetscCall(ISDestroy(&is[i])); 2887 2888 k = 0; 2889 for (j=0; j<ov; j++) { /* for each overlap */ 2890 n = isz; 2891 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2892 row = nidx[k]; 2893 start = ai[row]; 2894 end = ai[row+1]; 2895 for (l = start; l<end; l++) { 2896 val = aj[l]; 2897 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2898 } 2899 } 2900 } 2901 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i))); 2902 } 2903 PetscCall(PetscBTDestroy(&table)); 2904 PetscCall(PetscFree(nidx)); 2905 PetscFunctionReturn(0); 2906 } 2907 2908 /* -------------------------------------------------------------- */ 2909 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2910 { 2911 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2912 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2913 const PetscInt *row,*col; 2914 PetscInt *cnew,j,*lens; 2915 IS icolp,irowp; 2916 PetscInt *cwork = NULL; 2917 PetscScalar *vwork = NULL; 2918 2919 PetscFunctionBegin; 2920 PetscCall(ISInvertPermutation(rowp,PETSC_DECIDE,&irowp)); 2921 PetscCall(ISGetIndices(irowp,&row)); 2922 PetscCall(ISInvertPermutation(colp,PETSC_DECIDE,&icolp)); 2923 PetscCall(ISGetIndices(icolp,&col)); 2924 2925 /* determine lengths of permuted rows */ 2926 PetscCall(PetscMalloc1(m+1,&lens)); 2927 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2928 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 2929 PetscCall(MatSetSizes(*B,m,n,m,n)); 2930 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 2931 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2932 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens)); 2933 PetscCall(PetscFree(lens)); 2934 2935 PetscCall(PetscMalloc1(n,&cnew)); 2936 for (i=0; i<m; i++) { 2937 PetscCall(MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2938 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2939 PetscCall(MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES)); 2940 PetscCall(MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2941 } 2942 PetscCall(PetscFree(cnew)); 2943 2944 (*B)->assembled = PETSC_FALSE; 2945 2946 #if defined(PETSC_HAVE_DEVICE) 2947 PetscCall(MatBindToCPU(*B,A->boundtocpu)); 2948 #endif 2949 PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 2950 PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY)); 2951 PetscCall(ISRestoreIndices(irowp,&row)); 2952 PetscCall(ISRestoreIndices(icolp,&col)); 2953 PetscCall(ISDestroy(&irowp)); 2954 PetscCall(ISDestroy(&icolp)); 2955 if (rowp == colp) { 2956 PetscCall(MatPropagateSymmetryOptions(A,*B)); 2957 } 2958 PetscFunctionReturn(0); 2959 } 2960 2961 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2962 { 2963 PetscFunctionBegin; 2964 /* If the two matrices have the same copy implementation, use fast copy. */ 2965 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2966 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2967 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2968 const PetscScalar *aa; 2969 2970 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2971 PetscCheck(a->i[A->rmap->n] == b->i[B->rmap->n],PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different %" PetscInt_FMT " != %" PetscInt_FMT,a->i[A->rmap->n],b->i[B->rmap->n]); 2972 PetscCall(PetscArraycpy(b->a,aa,a->i[A->rmap->n])); 2973 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 2974 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2975 } else { 2976 PetscCall(MatCopy_Basic(A,B,str)); 2977 } 2978 PetscFunctionReturn(0); 2979 } 2980 2981 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2982 { 2983 PetscFunctionBegin; 2984 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,NULL)); 2985 PetscFunctionReturn(0); 2986 } 2987 2988 PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2989 { 2990 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2991 2992 PetscFunctionBegin; 2993 *array = a->a; 2994 PetscFunctionReturn(0); 2995 } 2996 2997 PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2998 { 2999 PetscFunctionBegin; 3000 *array = NULL; 3001 PetscFunctionReturn(0); 3002 } 3003 3004 /* 3005 Computes the number of nonzeros per row needed for preallocation when X and Y 3006 have different nonzero structure. 3007 */ 3008 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 3009 { 3010 PetscInt i,j,k,nzx,nzy; 3011 3012 PetscFunctionBegin; 3013 /* Set the number of nonzeros in the new matrix */ 3014 for (i=0; i<m; i++) { 3015 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 3016 nzx = xi[i+1] - xi[i]; 3017 nzy = yi[i+1] - yi[i]; 3018 nnz[i] = 0; 3019 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 3020 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 3021 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 3022 nnz[i]++; 3023 } 3024 for (; k<nzy; k++) nnz[i]++; 3025 } 3026 PetscFunctionReturn(0); 3027 } 3028 3029 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 3030 { 3031 PetscInt m = Y->rmap->N; 3032 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 3033 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 3034 3035 PetscFunctionBegin; 3036 /* Set the number of nonzeros in the new matrix */ 3037 PetscCall(MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz)); 3038 PetscFunctionReturn(0); 3039 } 3040 3041 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 3042 { 3043 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 3044 3045 PetscFunctionBegin; 3046 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 3047 PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE; 3048 if (e) { 3049 PetscCall(PetscArraycmp(x->i,y->i,Y->rmap->n+1,&e)); 3050 if (e) { 3051 PetscCall(PetscArraycmp(x->j,y->j,y->nz,&e)); 3052 if (e) str = SAME_NONZERO_PATTERN; 3053 } 3054 } 3055 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatStructure is not SAME_NONZERO_PATTERN"); 3056 } 3057 if (str == SAME_NONZERO_PATTERN) { 3058 const PetscScalar *xa; 3059 PetscScalar *ya,alpha = a; 3060 PetscBLASInt one = 1,bnz; 3061 3062 PetscCall(PetscBLASIntCast(x->nz,&bnz)); 3063 PetscCall(MatSeqAIJGetArray(Y,&ya)); 3064 PetscCall(MatSeqAIJGetArrayRead(X,&xa)); 3065 PetscCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,xa,&one,ya,&one)); 3066 PetscCall(MatSeqAIJRestoreArrayRead(X,&xa)); 3067 PetscCall(MatSeqAIJRestoreArray(Y,&ya)); 3068 PetscCall(PetscLogFlops(2.0*bnz)); 3069 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 3070 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3071 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 3072 PetscCall(MatAXPY_Basic(Y,a,X,str)); 3073 } else { 3074 Mat B; 3075 PetscInt *nnz; 3076 PetscCall(PetscMalloc1(Y->rmap->N,&nnz)); 3077 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),&B)); 3078 PetscCall(PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name)); 3079 PetscCall(MatSetLayouts(B,Y->rmap,Y->cmap)); 3080 PetscCall(MatSetType(B,((PetscObject)Y)->type_name)); 3081 PetscCall(MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz)); 3082 PetscCall(MatSeqAIJSetPreallocation(B,0,nnz)); 3083 PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 3084 PetscCall(MatHeaderMerge(Y,&B)); 3085 PetscCall(MatSeqAIJCheckInode(Y)); 3086 PetscCall(PetscFree(nnz)); 3087 } 3088 PetscFunctionReturn(0); 3089 } 3090 3091 PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 3092 { 3093 #if defined(PETSC_USE_COMPLEX) 3094 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3095 PetscInt i,nz; 3096 PetscScalar *a; 3097 3098 PetscFunctionBegin; 3099 nz = aij->nz; 3100 PetscCall(MatSeqAIJGetArray(mat,&a)); 3101 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 3102 PetscCall(MatSeqAIJRestoreArray(mat,&a)); 3103 #else 3104 PetscFunctionBegin; 3105 #endif 3106 PetscFunctionReturn(0); 3107 } 3108 3109 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3110 { 3111 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3112 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3113 PetscReal atmp; 3114 PetscScalar *x; 3115 const MatScalar *aa,*av; 3116 3117 PetscFunctionBegin; 3118 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3119 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3120 aa = av; 3121 ai = a->i; 3122 aj = a->j; 3123 3124 PetscCall(VecSet(v,0.0)); 3125 PetscCall(VecGetArrayWrite(v,&x)); 3126 PetscCall(VecGetLocalSize(v,&n)); 3127 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3128 for (i=0; i<m; i++) { 3129 ncols = ai[1] - ai[0]; ai++; 3130 for (j=0; j<ncols; j++) { 3131 atmp = PetscAbsScalar(*aa); 3132 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3133 aa++; aj++; 3134 } 3135 } 3136 PetscCall(VecRestoreArrayWrite(v,&x)); 3137 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3138 PetscFunctionReturn(0); 3139 } 3140 3141 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3142 { 3143 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3144 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3145 PetscScalar *x; 3146 const MatScalar *aa,*av; 3147 3148 PetscFunctionBegin; 3149 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3150 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3151 aa = av; 3152 ai = a->i; 3153 aj = a->j; 3154 3155 PetscCall(VecSet(v,0.0)); 3156 PetscCall(VecGetArrayWrite(v,&x)); 3157 PetscCall(VecGetLocalSize(v,&n)); 3158 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3159 for (i=0; i<m; i++) { 3160 ncols = ai[1] - ai[0]; ai++; 3161 if (ncols == A->cmap->n) { /* row is dense */ 3162 x[i] = *aa; if (idx) idx[i] = 0; 3163 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3164 x[i] = 0.0; 3165 if (idx) { 3166 for (j=0; j<ncols; j++) { /* find first implicit 0.0 in the row */ 3167 if (aj[j] > j) { 3168 idx[i] = j; 3169 break; 3170 } 3171 } 3172 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3173 if (j==ncols && j < A->cmap->n) idx[i] = j; 3174 } 3175 } 3176 for (j=0; j<ncols; j++) { 3177 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3178 aa++; aj++; 3179 } 3180 } 3181 PetscCall(VecRestoreArrayWrite(v,&x)); 3182 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3183 PetscFunctionReturn(0); 3184 } 3185 3186 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3187 { 3188 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3189 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3190 PetscScalar *x; 3191 const MatScalar *aa,*av; 3192 3193 PetscFunctionBegin; 3194 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3195 aa = av; 3196 ai = a->i; 3197 aj = a->j; 3198 3199 PetscCall(VecSet(v,0.0)); 3200 PetscCall(VecGetArrayWrite(v,&x)); 3201 PetscCall(VecGetLocalSize(v,&n)); 3202 PetscCheck(n == m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %" PetscInt_FMT " vs. %" PetscInt_FMT " rows", m, n); 3203 for (i=0; i<m; i++) { 3204 ncols = ai[1] - ai[0]; ai++; 3205 if (ncols == A->cmap->n) { /* row is dense */ 3206 x[i] = *aa; if (idx) idx[i] = 0; 3207 } else { /* row is sparse so already KNOW minimum is 0.0 or higher */ 3208 x[i] = 0.0; 3209 if (idx) { /* find first implicit 0.0 in the row */ 3210 for (j=0; j<ncols; j++) { 3211 if (aj[j] > j) { 3212 idx[i] = j; 3213 break; 3214 } 3215 } 3216 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3217 if (j==ncols && j < A->cmap->n) idx[i] = j; 3218 } 3219 } 3220 for (j=0; j<ncols; j++) { 3221 if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3222 aa++; aj++; 3223 } 3224 } 3225 PetscCall(VecRestoreArrayWrite(v,&x)); 3226 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3227 PetscFunctionReturn(0); 3228 } 3229 3230 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3231 { 3232 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3233 PetscInt i,j,m = A->rmap->n,ncols,n; 3234 const PetscInt *ai,*aj; 3235 PetscScalar *x; 3236 const MatScalar *aa,*av; 3237 3238 PetscFunctionBegin; 3239 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3240 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3241 aa = av; 3242 ai = a->i; 3243 aj = a->j; 3244 3245 PetscCall(VecSet(v,0.0)); 3246 PetscCall(VecGetArrayWrite(v,&x)); 3247 PetscCall(VecGetLocalSize(v,&n)); 3248 PetscCheck(n == m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3249 for (i=0; i<m; i++) { 3250 ncols = ai[1] - ai[0]; ai++; 3251 if (ncols == A->cmap->n) { /* row is dense */ 3252 x[i] = *aa; if (idx) idx[i] = 0; 3253 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3254 x[i] = 0.0; 3255 if (idx) { /* find first implicit 0.0 in the row */ 3256 for (j=0; j<ncols; j++) { 3257 if (aj[j] > j) { 3258 idx[i] = j; 3259 break; 3260 } 3261 } 3262 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3263 if (j==ncols && j < A->cmap->n) idx[i] = j; 3264 } 3265 } 3266 for (j=0; j<ncols; j++) { 3267 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3268 aa++; aj++; 3269 } 3270 } 3271 PetscCall(VecRestoreArrayWrite(v,&x)); 3272 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3273 PetscFunctionReturn(0); 3274 } 3275 3276 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3277 { 3278 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3279 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3280 MatScalar *diag,work[25],*v_work; 3281 const PetscReal shift = 0.0; 3282 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3283 3284 PetscFunctionBegin; 3285 allowzeropivot = PetscNot(A->erroriffailure); 3286 if (a->ibdiagvalid) { 3287 if (values) *values = a->ibdiag; 3288 PetscFunctionReturn(0); 3289 } 3290 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 3291 if (!a->ibdiag) { 3292 PetscCall(PetscMalloc1(bs2*mbs,&a->ibdiag)); 3293 PetscCall(PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar))); 3294 } 3295 diag = a->ibdiag; 3296 if (values) *values = a->ibdiag; 3297 /* factor and invert each block */ 3298 switch (bs) { 3299 case 1: 3300 for (i=0; i<mbs; i++) { 3301 PetscCall(MatGetValues(A,1,&i,1,&i,diag+i)); 3302 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3303 if (allowzeropivot) { 3304 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3305 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3306 A->factorerror_zeropivot_row = i; 3307 PetscCall(PetscInfo(A,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON)); 3308 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON); 3309 } 3310 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3311 } 3312 break; 3313 case 2: 3314 for (i=0; i<mbs; i++) { 3315 ij[0] = 2*i; ij[1] = 2*i + 1; 3316 PetscCall(MatGetValues(A,2,ij,2,ij,diag)); 3317 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 3318 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3319 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 3320 diag += 4; 3321 } 3322 break; 3323 case 3: 3324 for (i=0; i<mbs; i++) { 3325 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3326 PetscCall(MatGetValues(A,3,ij,3,ij,diag)); 3327 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 3328 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3329 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 3330 diag += 9; 3331 } 3332 break; 3333 case 4: 3334 for (i=0; i<mbs; i++) { 3335 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3336 PetscCall(MatGetValues(A,4,ij,4,ij,diag)); 3337 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 3338 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3339 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 3340 diag += 16; 3341 } 3342 break; 3343 case 5: 3344 for (i=0; i<mbs; i++) { 3345 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3346 PetscCall(MatGetValues(A,5,ij,5,ij,diag)); 3347 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 3348 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3349 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 3350 diag += 25; 3351 } 3352 break; 3353 case 6: 3354 for (i=0; i<mbs; i++) { 3355 ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5; 3356 PetscCall(MatGetValues(A,6,ij,6,ij,diag)); 3357 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 3358 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3359 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 3360 diag += 36; 3361 } 3362 break; 3363 case 7: 3364 for (i=0; i<mbs; i++) { 3365 ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6; 3366 PetscCall(MatGetValues(A,7,ij,7,ij,diag)); 3367 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 3368 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3369 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 3370 diag += 49; 3371 } 3372 break; 3373 default: 3374 PetscCall(PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ)); 3375 for (i=0; i<mbs; i++) { 3376 for (j=0; j<bs; j++) { 3377 IJ[j] = bs*i + j; 3378 } 3379 PetscCall(MatGetValues(A,bs,IJ,bs,IJ,diag)); 3380 PetscCall(PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 3381 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3382 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bs)); 3383 diag += bs2; 3384 } 3385 PetscCall(PetscFree3(v_work,v_pivots,IJ)); 3386 } 3387 a->ibdiagvalid = PETSC_TRUE; 3388 PetscFunctionReturn(0); 3389 } 3390 3391 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3392 { 3393 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3394 PetscScalar a,*aa; 3395 PetscInt m,n,i,j,col; 3396 3397 PetscFunctionBegin; 3398 if (!x->assembled) { 3399 PetscCall(MatGetSize(x,&m,&n)); 3400 for (i=0; i<m; i++) { 3401 for (j=0; j<aij->imax[i]; j++) { 3402 PetscCall(PetscRandomGetValue(rctx,&a)); 3403 col = (PetscInt)(n*PetscRealPart(a)); 3404 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3405 } 3406 } 3407 } else { 3408 PetscCall(MatSeqAIJGetArrayWrite(x,&aa)); 3409 for (i=0; i<aij->nz; i++) PetscCall(PetscRandomGetValue(rctx,aa+i)); 3410 PetscCall(MatSeqAIJRestoreArrayWrite(x,&aa)); 3411 } 3412 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3413 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3414 PetscFunctionReturn(0); 3415 } 3416 3417 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3418 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx) 3419 { 3420 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3421 PetscScalar a; 3422 PetscInt m,n,i,j,col,nskip; 3423 3424 PetscFunctionBegin; 3425 nskip = high - low; 3426 PetscCall(MatGetSize(x,&m,&n)); 3427 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3428 for (i=0; i<m; i++) { 3429 for (j=0; j<aij->imax[i]; j++) { 3430 PetscCall(PetscRandomGetValue(rctx,&a)); 3431 col = (PetscInt)(n*PetscRealPart(a)); 3432 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3433 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3434 } 3435 } 3436 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3437 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3438 PetscFunctionReturn(0); 3439 } 3440 3441 /* -------------------------------------------------------------------*/ 3442 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3443 MatGetRow_SeqAIJ, 3444 MatRestoreRow_SeqAIJ, 3445 MatMult_SeqAIJ, 3446 /* 4*/ MatMultAdd_SeqAIJ, 3447 MatMultTranspose_SeqAIJ, 3448 MatMultTransposeAdd_SeqAIJ, 3449 NULL, 3450 NULL, 3451 NULL, 3452 /* 10*/ NULL, 3453 MatLUFactor_SeqAIJ, 3454 NULL, 3455 MatSOR_SeqAIJ, 3456 MatTranspose_SeqAIJ, 3457 /*1 5*/ MatGetInfo_SeqAIJ, 3458 MatEqual_SeqAIJ, 3459 MatGetDiagonal_SeqAIJ, 3460 MatDiagonalScale_SeqAIJ, 3461 MatNorm_SeqAIJ, 3462 /* 20*/ NULL, 3463 MatAssemblyEnd_SeqAIJ, 3464 MatSetOption_SeqAIJ, 3465 MatZeroEntries_SeqAIJ, 3466 /* 24*/ MatZeroRows_SeqAIJ, 3467 NULL, 3468 NULL, 3469 NULL, 3470 NULL, 3471 /* 29*/ MatSetUp_SeqAIJ, 3472 NULL, 3473 NULL, 3474 NULL, 3475 NULL, 3476 /* 34*/ MatDuplicate_SeqAIJ, 3477 NULL, 3478 NULL, 3479 MatILUFactor_SeqAIJ, 3480 NULL, 3481 /* 39*/ MatAXPY_SeqAIJ, 3482 MatCreateSubMatrices_SeqAIJ, 3483 MatIncreaseOverlap_SeqAIJ, 3484 MatGetValues_SeqAIJ, 3485 MatCopy_SeqAIJ, 3486 /* 44*/ MatGetRowMax_SeqAIJ, 3487 MatScale_SeqAIJ, 3488 MatShift_SeqAIJ, 3489 MatDiagonalSet_SeqAIJ, 3490 MatZeroRowsColumns_SeqAIJ, 3491 /* 49*/ MatSetRandom_SeqAIJ, 3492 MatGetRowIJ_SeqAIJ, 3493 MatRestoreRowIJ_SeqAIJ, 3494 MatGetColumnIJ_SeqAIJ, 3495 MatRestoreColumnIJ_SeqAIJ, 3496 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3497 NULL, 3498 NULL, 3499 MatPermute_SeqAIJ, 3500 NULL, 3501 /* 59*/ NULL, 3502 MatDestroy_SeqAIJ, 3503 MatView_SeqAIJ, 3504 NULL, 3505 NULL, 3506 /* 64*/ NULL, 3507 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3508 NULL, 3509 NULL, 3510 NULL, 3511 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3512 MatGetRowMinAbs_SeqAIJ, 3513 NULL, 3514 NULL, 3515 NULL, 3516 /* 74*/ NULL, 3517 MatFDColoringApply_AIJ, 3518 NULL, 3519 NULL, 3520 NULL, 3521 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3522 NULL, 3523 NULL, 3524 NULL, 3525 MatLoad_SeqAIJ, 3526 /* 84*/ MatIsSymmetric_SeqAIJ, 3527 MatIsHermitian_SeqAIJ, 3528 NULL, 3529 NULL, 3530 NULL, 3531 /* 89*/ NULL, 3532 NULL, 3533 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3534 NULL, 3535 NULL, 3536 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3537 NULL, 3538 NULL, 3539 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3540 NULL, 3541 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3542 NULL, 3543 NULL, 3544 MatConjugate_SeqAIJ, 3545 NULL, 3546 /*104*/ MatSetValuesRow_SeqAIJ, 3547 MatRealPart_SeqAIJ, 3548 MatImaginaryPart_SeqAIJ, 3549 NULL, 3550 NULL, 3551 /*109*/ MatMatSolve_SeqAIJ, 3552 NULL, 3553 MatGetRowMin_SeqAIJ, 3554 NULL, 3555 MatMissingDiagonal_SeqAIJ, 3556 /*114*/ NULL, 3557 NULL, 3558 NULL, 3559 NULL, 3560 NULL, 3561 /*119*/ NULL, 3562 NULL, 3563 NULL, 3564 NULL, 3565 MatGetMultiProcBlock_SeqAIJ, 3566 /*124*/ MatFindNonzeroRows_SeqAIJ, 3567 MatGetColumnReductions_SeqAIJ, 3568 MatInvertBlockDiagonal_SeqAIJ, 3569 MatInvertVariableBlockDiagonal_SeqAIJ, 3570 NULL, 3571 /*129*/ NULL, 3572 NULL, 3573 NULL, 3574 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3575 MatTransposeColoringCreate_SeqAIJ, 3576 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3577 MatTransColoringApplyDenToSp_SeqAIJ, 3578 NULL, 3579 NULL, 3580 MatRARtNumeric_SeqAIJ_SeqAIJ, 3581 /*139*/NULL, 3582 NULL, 3583 NULL, 3584 MatFDColoringSetUp_SeqXAIJ, 3585 MatFindOffBlockDiagonalEntries_SeqAIJ, 3586 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3587 /*145*/MatDestroySubMatrices_SeqAIJ, 3588 NULL, 3589 NULL, 3590 MatCreateGraph_Simple_AIJ, 3591 MatFilter_AIJ, 3592 /*150*/MatTransposeSymbolic_SeqAIJ 3593 }; 3594 3595 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3596 { 3597 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3598 PetscInt i,nz,n; 3599 3600 PetscFunctionBegin; 3601 nz = aij->maxnz; 3602 n = mat->rmap->n; 3603 for (i=0; i<nz; i++) { 3604 aij->j[i] = indices[i]; 3605 } 3606 aij->nz = nz; 3607 for (i=0; i<n; i++) { 3608 aij->ilen[i] = aij->imax[i]; 3609 } 3610 PetscFunctionReturn(0); 3611 } 3612 3613 /* 3614 * Given a sparse matrix with global column indices, compact it by using a local column space. 3615 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3616 */ 3617 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3618 { 3619 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3620 PetscTable gid1_lid1; 3621 PetscTablePosition tpos; 3622 PetscInt gid,lid,i,ec,nz = aij->nz; 3623 PetscInt *garray,*jj = aij->j; 3624 3625 PetscFunctionBegin; 3626 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3627 PetscValidPointer(mapping,2); 3628 /* use a table */ 3629 PetscCall(PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1)); 3630 ec = 0; 3631 for (i=0; i<nz; i++) { 3632 PetscInt data,gid1 = jj[i] + 1; 3633 PetscCall(PetscTableFind(gid1_lid1,gid1,&data)); 3634 if (!data) { 3635 /* one based table */ 3636 PetscCall(PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES)); 3637 } 3638 } 3639 /* form array of columns we need */ 3640 PetscCall(PetscMalloc1(ec,&garray)); 3641 PetscCall(PetscTableGetHeadPosition(gid1_lid1,&tpos)); 3642 while (tpos) { 3643 PetscCall(PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid)); 3644 gid--; 3645 lid--; 3646 garray[lid] = gid; 3647 } 3648 PetscCall(PetscSortInt(ec,garray)); /* sort, and rebuild */ 3649 PetscCall(PetscTableRemoveAll(gid1_lid1)); 3650 for (i=0; i<ec; i++) { 3651 PetscCall(PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES)); 3652 } 3653 /* compact out the extra columns in B */ 3654 for (i=0; i<nz; i++) { 3655 PetscInt gid1 = jj[i] + 1; 3656 PetscCall(PetscTableFind(gid1_lid1,gid1,&lid)); 3657 lid--; 3658 jj[i] = lid; 3659 } 3660 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3661 PetscCall(PetscTableDestroy(&gid1_lid1)); 3662 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap)); 3663 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping)); 3664 PetscCall(ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH)); 3665 PetscFunctionReturn(0); 3666 } 3667 3668 /*@ 3669 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3670 in the matrix. 3671 3672 Input Parameters: 3673 + mat - the SeqAIJ matrix 3674 - indices - the column indices 3675 3676 Level: advanced 3677 3678 Notes: 3679 This can be called if you have precomputed the nonzero structure of the 3680 matrix and want to provide it to the matrix object to improve the performance 3681 of the MatSetValues() operation. 3682 3683 You MUST have set the correct numbers of nonzeros per row in the call to 3684 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3685 3686 MUST be called before any calls to MatSetValues(); 3687 3688 The indices should start with zero, not one. 3689 3690 @*/ 3691 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3692 { 3693 PetscFunctionBegin; 3694 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3695 PetscValidIntPointer(indices,2); 3696 PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices)); 3697 PetscFunctionReturn(0); 3698 } 3699 3700 /* ----------------------------------------------------------------------------------------*/ 3701 3702 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3703 { 3704 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3705 size_t nz = aij->i[mat->rmap->n]; 3706 3707 PetscFunctionBegin; 3708 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3709 3710 /* allocate space for values if not already there */ 3711 if (!aij->saved_values) { 3712 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 3713 PetscCall(PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar))); 3714 } 3715 3716 /* copy values over */ 3717 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 3718 PetscFunctionReturn(0); 3719 } 3720 3721 /*@ 3722 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3723 example, reuse of the linear part of a Jacobian, while recomputing the 3724 nonlinear portion. 3725 3726 Collect on Mat 3727 3728 Input Parameters: 3729 . mat - the matrix (currently only AIJ matrices support this option) 3730 3731 Level: advanced 3732 3733 Common Usage, with SNESSolve(): 3734 $ Create Jacobian matrix 3735 $ Set linear terms into matrix 3736 $ Apply boundary conditions to matrix, at this time matrix must have 3737 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3738 $ boundary conditions again will not change the nonzero structure 3739 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3740 $ ierr = MatStoreValues(mat); 3741 $ Call SNESSetJacobian() with matrix 3742 $ In your Jacobian routine 3743 $ ierr = MatRetrieveValues(mat); 3744 $ Set nonlinear terms in matrix 3745 3746 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3747 $ // build linear portion of Jacobian 3748 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3749 $ ierr = MatStoreValues(mat); 3750 $ loop over nonlinear iterations 3751 $ ierr = MatRetrieveValues(mat); 3752 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3753 $ // call MatAssemblyBegin/End() on matrix 3754 $ Solve linear system with Jacobian 3755 $ endloop 3756 3757 Notes: 3758 Matrix must already be assemblied before calling this routine 3759 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3760 calling this routine. 3761 3762 When this is called multiple times it overwrites the previous set of stored values 3763 and does not allocated additional space. 3764 3765 .seealso: `MatRetrieveValues()` 3766 3767 @*/ 3768 PetscErrorCode MatStoreValues(Mat mat) 3769 { 3770 PetscFunctionBegin; 3771 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3772 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3773 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3774 PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat)); 3775 PetscFunctionReturn(0); 3776 } 3777 3778 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3779 { 3780 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3781 PetscInt nz = aij->i[mat->rmap->n]; 3782 3783 PetscFunctionBegin; 3784 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3785 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3786 /* copy values over */ 3787 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 3788 PetscFunctionReturn(0); 3789 } 3790 3791 /*@ 3792 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3793 example, reuse of the linear part of a Jacobian, while recomputing the 3794 nonlinear portion. 3795 3796 Collect on Mat 3797 3798 Input Parameters: 3799 . mat - the matrix (currently only AIJ matrices support this option) 3800 3801 Level: advanced 3802 3803 .seealso: `MatStoreValues()` 3804 3805 @*/ 3806 PetscErrorCode MatRetrieveValues(Mat mat) 3807 { 3808 PetscFunctionBegin; 3809 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3810 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3811 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3812 PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat)); 3813 PetscFunctionReturn(0); 3814 } 3815 3816 /* --------------------------------------------------------------------------------*/ 3817 /*@C 3818 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3819 (the default parallel PETSc format). For good matrix assembly performance 3820 the user should preallocate the matrix storage by setting the parameter nz 3821 (or the array nnz). By setting these parameters accurately, performance 3822 during matrix assembly can be increased by more than a factor of 50. 3823 3824 Collective 3825 3826 Input Parameters: 3827 + comm - MPI communicator, set to PETSC_COMM_SELF 3828 . m - number of rows 3829 . n - number of columns 3830 . nz - number of nonzeros per row (same for all rows) 3831 - nnz - array containing the number of nonzeros in the various rows 3832 (possibly different for each row) or NULL 3833 3834 Output Parameter: 3835 . A - the matrix 3836 3837 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3838 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3839 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3840 3841 Notes: 3842 If nnz is given then nz is ignored 3843 3844 The AIJ format (also called the Yale sparse matrix format or 3845 compressed row storage), is fully compatible with standard Fortran 77 3846 storage. That is, the stored row and column indices can begin at 3847 either one (as in Fortran) or zero. See the users' manual for details. 3848 3849 Specify the preallocated storage with either nz or nnz (not both). 3850 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3851 allocation. For large problems you MUST preallocate memory or you 3852 will get TERRIBLE performance, see the users' manual chapter on matrices. 3853 3854 By default, this format uses inodes (identical nodes) when possible, to 3855 improve numerical efficiency of matrix-vector products and solves. We 3856 search for consecutive rows with the same nonzero structure, thereby 3857 reusing matrix information to achieve increased efficiency. 3858 3859 Options Database Keys: 3860 + -mat_no_inode - Do not use inodes 3861 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3862 3863 Level: intermediate 3864 3865 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()` 3866 3867 @*/ 3868 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3869 { 3870 PetscFunctionBegin; 3871 PetscCall(MatCreate(comm,A)); 3872 PetscCall(MatSetSizes(*A,m,n,m,n)); 3873 PetscCall(MatSetType(*A,MATSEQAIJ)); 3874 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz)); 3875 PetscFunctionReturn(0); 3876 } 3877 3878 /*@C 3879 MatSeqAIJSetPreallocation - For good matrix assembly performance 3880 the user should preallocate the matrix storage by setting the parameter nz 3881 (or the array nnz). By setting these parameters accurately, performance 3882 during matrix assembly can be increased by more than a factor of 50. 3883 3884 Collective 3885 3886 Input Parameters: 3887 + B - The matrix 3888 . nz - number of nonzeros per row (same for all rows) 3889 - nnz - array containing the number of nonzeros in the various rows 3890 (possibly different for each row) or NULL 3891 3892 Notes: 3893 If nnz is given then nz is ignored 3894 3895 The AIJ format (also called the Yale sparse matrix format or 3896 compressed row storage), is fully compatible with standard Fortran 77 3897 storage. That is, the stored row and column indices can begin at 3898 either one (as in Fortran) or zero. See the users' manual for details. 3899 3900 Specify the preallocated storage with either nz or nnz (not both). 3901 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3902 allocation. For large problems you MUST preallocate memory or you 3903 will get TERRIBLE performance, see the users' manual chapter on matrices. 3904 3905 You can call MatGetInfo() to get information on how effective the preallocation was; 3906 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3907 You can also run with the option -info and look for messages with the string 3908 malloc in them to see if additional memory allocation was needed. 3909 3910 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3911 entries or columns indices 3912 3913 By default, this format uses inodes (identical nodes) when possible, to 3914 improve numerical efficiency of matrix-vector products and solves. We 3915 search for consecutive rows with the same nonzero structure, thereby 3916 reusing matrix information to achieve increased efficiency. 3917 3918 Options Database Keys: 3919 + -mat_no_inode - Do not use inodes 3920 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3921 3922 Level: intermediate 3923 3924 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatGetInfo()`, 3925 `MatSeqAIJSetTotalPreallocation()` 3926 3927 @*/ 3928 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3929 { 3930 PetscFunctionBegin; 3931 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3932 PetscValidType(B,1); 3933 PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz)); 3934 PetscFunctionReturn(0); 3935 } 3936 3937 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3938 { 3939 Mat_SeqAIJ *b; 3940 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3941 PetscInt i; 3942 3943 PetscFunctionBegin; 3944 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3945 if (nz == MAT_SKIP_ALLOCATION) { 3946 skipallocation = PETSC_TRUE; 3947 nz = 0; 3948 } 3949 PetscCall(PetscLayoutSetUp(B->rmap)); 3950 PetscCall(PetscLayoutSetUp(B->cmap)); 3951 3952 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3953 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 3954 if (PetscUnlikelyDebug(nnz)) { 3955 for (i=0; i<B->rmap->n; i++) { 3956 PetscCheck(nnz[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT,i,nnz[i]); 3957 PetscCheck(nnz[i] <= B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT,i,nnz[i],B->cmap->n); 3958 } 3959 } 3960 3961 B->preallocated = PETSC_TRUE; 3962 3963 b = (Mat_SeqAIJ*)B->data; 3964 3965 if (!skipallocation) { 3966 if (!b->imax) { 3967 PetscCall(PetscMalloc1(B->rmap->n,&b->imax)); 3968 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 3969 } 3970 if (!b->ilen) { 3971 /* b->ilen will count nonzeros in each row so far. */ 3972 PetscCall(PetscCalloc1(B->rmap->n,&b->ilen)); 3973 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 3974 } else { 3975 PetscCall(PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt))); 3976 } 3977 if (!b->ipre) { 3978 PetscCall(PetscMalloc1(B->rmap->n,&b->ipre)); 3979 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 3980 } 3981 if (!nnz) { 3982 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3983 else if (nz < 0) nz = 1; 3984 nz = PetscMin(nz,B->cmap->n); 3985 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3986 nz = nz*B->rmap->n; 3987 } else { 3988 PetscInt64 nz64 = 0; 3989 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 3990 PetscCall(PetscIntCast(nz64,&nz)); 3991 } 3992 3993 /* allocate the matrix space */ 3994 /* FIXME: should B's old memory be unlogged? */ 3995 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 3996 if (B->structure_only) { 3997 PetscCall(PetscMalloc1(nz,&b->j)); 3998 PetscCall(PetscMalloc1(B->rmap->n+1,&b->i)); 3999 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt))); 4000 } else { 4001 PetscCall(PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i)); 4002 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)))); 4003 } 4004 b->i[0] = 0; 4005 for (i=1; i<B->rmap->n+1; i++) { 4006 b->i[i] = b->i[i-1] + b->imax[i-1]; 4007 } 4008 if (B->structure_only) { 4009 b->singlemalloc = PETSC_FALSE; 4010 b->free_a = PETSC_FALSE; 4011 } else { 4012 b->singlemalloc = PETSC_TRUE; 4013 b->free_a = PETSC_TRUE; 4014 } 4015 b->free_ij = PETSC_TRUE; 4016 } else { 4017 b->free_a = PETSC_FALSE; 4018 b->free_ij = PETSC_FALSE; 4019 } 4020 4021 if (b->ipre && nnz != b->ipre && b->imax) { 4022 /* reserve user-requested sparsity */ 4023 PetscCall(PetscArraycpy(b->ipre,b->imax,B->rmap->n)); 4024 } 4025 4026 b->nz = 0; 4027 b->maxnz = nz; 4028 B->info.nz_unneeded = (double)b->maxnz; 4029 if (realalloc) PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 4030 B->was_assembled = PETSC_FALSE; 4031 B->assembled = PETSC_FALSE; 4032 /* We simply deem preallocation has changed nonzero state. Updating the state 4033 will give clients (like AIJKokkos) a chance to know something has happened. 4034 */ 4035 B->nonzerostate++; 4036 PetscFunctionReturn(0); 4037 } 4038 4039 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4040 { 4041 Mat_SeqAIJ *a; 4042 PetscInt i; 4043 4044 PetscFunctionBegin; 4045 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4046 4047 /* Check local size. If zero, then return */ 4048 if (!A->rmap->n) PetscFunctionReturn(0); 4049 4050 a = (Mat_SeqAIJ*)A->data; 4051 /* if no saved info, we error out */ 4052 PetscCheck(a->ipre,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"No saved preallocation info "); 4053 4054 PetscCheck(a->i && a->j && a->a && a->imax && a->ilen,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Memory info is incomplete, and can not reset preallocation "); 4055 4056 PetscCall(PetscArraycpy(a->imax,a->ipre,A->rmap->n)); 4057 PetscCall(PetscArrayzero(a->ilen,A->rmap->n)); 4058 a->i[0] = 0; 4059 for (i=1; i<A->rmap->n+1; i++) { 4060 a->i[i] = a->i[i-1] + a->imax[i-1]; 4061 } 4062 A->preallocated = PETSC_TRUE; 4063 a->nz = 0; 4064 a->maxnz = a->i[A->rmap->n]; 4065 A->info.nz_unneeded = (double)a->maxnz; 4066 A->was_assembled = PETSC_FALSE; 4067 A->assembled = PETSC_FALSE; 4068 PetscFunctionReturn(0); 4069 } 4070 4071 /*@ 4072 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 4073 4074 Input Parameters: 4075 + B - the matrix 4076 . i - the indices into j for the start of each row (starts with zero) 4077 . j - the column indices for each row (starts with zero) these must be sorted for each row 4078 - v - optional values in the matrix 4079 4080 Level: developer 4081 4082 Notes: 4083 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 4084 4085 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4086 structure will be the union of all the previous nonzero structures. 4087 4088 Developer Notes: 4089 An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and 4090 then just copies the v values directly with PetscMemcpy(). 4091 4092 This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them. 4093 4094 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatSeqAIJSetPreallocation()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MatResetPreallocation()` 4095 @*/ 4096 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 4097 { 4098 PetscFunctionBegin; 4099 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 4100 PetscValidType(B,1); 4101 PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v)); 4102 PetscFunctionReturn(0); 4103 } 4104 4105 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 4106 { 4107 PetscInt i; 4108 PetscInt m,n; 4109 PetscInt nz; 4110 PetscInt *nnz; 4111 4112 PetscFunctionBegin; 4113 PetscCheck(Ii[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4114 4115 PetscCall(PetscLayoutSetUp(B->rmap)); 4116 PetscCall(PetscLayoutSetUp(B->cmap)); 4117 4118 PetscCall(MatGetSize(B, &m, &n)); 4119 PetscCall(PetscMalloc1(m+1, &nnz)); 4120 for (i = 0; i < m; i++) { 4121 nz = Ii[i+1]- Ii[i]; 4122 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4123 nnz[i] = nz; 4124 } 4125 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4126 PetscCall(PetscFree(nnz)); 4127 4128 for (i = 0; i < m; i++) { 4129 PetscCall(MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES)); 4130 } 4131 4132 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 4133 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 4134 4135 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 4136 PetscFunctionReturn(0); 4137 } 4138 4139 /*@ 4140 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4141 4142 Input Parameters: 4143 + A - left-hand side matrix 4144 . B - right-hand side matrix 4145 - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4146 4147 Output Parameter: 4148 . C - Kronecker product of A and B 4149 4150 Level: intermediate 4151 4152 Notes: 4153 MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron(). 4154 4155 .seealso: `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATKAIJ`, `MatReuse` 4156 @*/ 4157 PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C) 4158 { 4159 PetscFunctionBegin; 4160 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4161 PetscValidType(A,1); 4162 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 4163 PetscValidType(B,2); 4164 PetscValidPointer(C,4); 4165 if (reuse == MAT_REUSE_MATRIX) { 4166 PetscValidHeaderSpecific(*C,MAT_CLASSID,4); 4167 PetscValidType(*C,4); 4168 } 4169 PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C)); 4170 PetscFunctionReturn(0); 4171 } 4172 4173 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C) 4174 { 4175 Mat newmat; 4176 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4177 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 4178 PetscScalar *v; 4179 const PetscScalar *aa,*ba; 4180 PetscInt *i,*j,m,n,p,q,nnz = 0,am = A->rmap->n,bm = B->rmap->n,an = A->cmap->n, bn = B->cmap->n; 4181 PetscBool flg; 4182 4183 PetscFunctionBegin; 4184 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4185 PetscCheck(A->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4186 PetscCheck(!B->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4187 PetscCheck(B->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4188 PetscCall(PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg)); 4189 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatType %s",((PetscObject)B)->type_name); 4190 PetscCheck(reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatReuse %d",(int)reuse); 4191 if (reuse == MAT_INITIAL_MATRIX) { 4192 PetscCall(PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j)); 4193 PetscCall(MatCreate(PETSC_COMM_SELF,&newmat)); 4194 PetscCall(MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn)); 4195 PetscCall(MatSetType(newmat,MATAIJ)); 4196 i[0] = 0; 4197 for (m = 0; m < am; ++m) { 4198 for (p = 0; p < bm; ++p) { 4199 i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]); 4200 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4201 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4202 j[nnz++] = a->j[n]*bn + b->j[q]; 4203 } 4204 } 4205 } 4206 } 4207 PetscCall(MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL)); 4208 *C = newmat; 4209 PetscCall(PetscFree2(i,j)); 4210 nnz = 0; 4211 } 4212 PetscCall(MatSeqAIJGetArray(*C,&v)); 4213 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4214 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 4215 for (m = 0; m < am; ++m) { 4216 for (p = 0; p < bm; ++p) { 4217 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4218 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4219 v[nnz++] = aa[n] * ba[q]; 4220 } 4221 } 4222 } 4223 } 4224 PetscCall(MatSeqAIJRestoreArray(*C,&v)); 4225 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 4226 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 4227 PetscFunctionReturn(0); 4228 } 4229 4230 #include <../src/mat/impls/dense/seq/dense.h> 4231 #include <petsc/private/kernels/petscaxpy.h> 4232 4233 /* 4234 Computes (B'*A')' since computing B*A directly is untenable 4235 4236 n p p 4237 [ ] [ ] [ ] 4238 m [ A ] * n [ B ] = m [ C ] 4239 [ ] [ ] [ ] 4240 4241 */ 4242 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4243 { 4244 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4245 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4246 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4247 PetscInt i,j,n,m,q,p; 4248 const PetscInt *ii,*idx; 4249 const PetscScalar *b,*a,*a_q; 4250 PetscScalar *c,*c_q; 4251 PetscInt clda = sub_c->lda; 4252 PetscInt alda = sub_a->lda; 4253 4254 PetscFunctionBegin; 4255 m = A->rmap->n; 4256 n = A->cmap->n; 4257 p = B->cmap->n; 4258 a = sub_a->v; 4259 b = sub_b->a; 4260 c = sub_c->v; 4261 if (clda == m) { 4262 PetscCall(PetscArrayzero(c,m*p)); 4263 } else { 4264 for (j=0;j<p;j++) 4265 for (i=0;i<m;i++) 4266 c[j*clda + i] = 0.0; 4267 } 4268 ii = sub_b->i; 4269 idx = sub_b->j; 4270 for (i=0; i<n; i++) { 4271 q = ii[i+1] - ii[i]; 4272 while (q-->0) { 4273 c_q = c + clda*(*idx); 4274 a_q = a + alda*i; 4275 PetscKernelAXPY(c_q,*b,a_q,m); 4276 idx++; 4277 b++; 4278 } 4279 } 4280 PetscFunctionReturn(0); 4281 } 4282 4283 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C) 4284 { 4285 PetscInt m=A->rmap->n,n=B->cmap->n; 4286 PetscBool cisdense; 4287 4288 PetscFunctionBegin; 4289 PetscCheck(A->cmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %" PetscInt_FMT " != B->rmap->n %" PetscInt_FMT,A->cmap->n,B->rmap->n); 4290 PetscCall(MatSetSizes(C,m,n,m,n)); 4291 PetscCall(MatSetBlockSizesFromMats(C,A,B)); 4292 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 4293 if (!cisdense) { 4294 PetscCall(MatSetType(C,MATDENSE)); 4295 } 4296 PetscCall(MatSetUp(C)); 4297 4298 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4299 PetscFunctionReturn(0); 4300 } 4301 4302 /* ----------------------------------------------------------------*/ 4303 /*MC 4304 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4305 based on compressed sparse row format. 4306 4307 Options Database Keys: 4308 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4309 4310 Level: beginner 4311 4312 Notes: 4313 MatSetValues() may be called for this matrix type with a NULL argument for the numerical values, 4314 in this case the values associated with the rows and columns one passes in are set to zero 4315 in the matrix 4316 4317 MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no 4318 space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored 4319 4320 Developer Notes: 4321 It would be nice if all matrix formats supported passing NULL in for the numerical values 4322 4323 .seealso: `MatCreateSeqAIJ()`, `MatSetFromOptions()`, `MatSetType()`, `MatCreate()`, `MatType`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4324 M*/ 4325 4326 /*MC 4327 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4328 4329 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4330 and MATMPIAIJ otherwise. As a result, for single process communicators, 4331 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4332 for communicators controlling multiple processes. It is recommended that you call both of 4333 the above preallocation routines for simplicity. 4334 4335 Options Database Keys: 4336 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4337 4338 Developer Notes: 4339 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4340 enough exist. 4341 4342 Level: beginner 4343 4344 .seealso: `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATMPIAIJ`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4345 M*/ 4346 4347 /*MC 4348 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4349 4350 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4351 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4352 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4353 for communicators controlling multiple processes. It is recommended that you call both of 4354 the above preallocation routines for simplicity. 4355 4356 Options Database Keys: 4357 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4358 4359 Level: beginner 4360 4361 .seealso: `MatCreateMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL` 4362 M*/ 4363 4364 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4365 #if defined(PETSC_HAVE_ELEMENTAL) 4366 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4367 #endif 4368 #if defined(PETSC_HAVE_SCALAPACK) 4369 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 4370 #endif 4371 #if defined(PETSC_HAVE_HYPRE) 4372 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4373 #endif 4374 4375 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4376 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4377 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4378 4379 /*@C 4380 MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored 4381 4382 Not Collective 4383 4384 Input Parameter: 4385 . mat - a MATSEQAIJ matrix 4386 4387 Output Parameter: 4388 . array - pointer to the data 4389 4390 Level: intermediate 4391 4392 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4393 @*/ 4394 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4395 { 4396 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4397 4398 PetscFunctionBegin; 4399 if (aij->ops->getarray) { 4400 PetscCall((*aij->ops->getarray)(A,array)); 4401 } else { 4402 *array = aij->a; 4403 } 4404 PetscFunctionReturn(0); 4405 } 4406 4407 /*@C 4408 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4409 4410 Not Collective 4411 4412 Input Parameters: 4413 + mat - a MATSEQAIJ matrix 4414 - array - pointer to the data 4415 4416 Level: intermediate 4417 4418 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayF90()` 4419 @*/ 4420 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4421 { 4422 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4423 4424 PetscFunctionBegin; 4425 if (aij->ops->restorearray) { 4426 PetscCall((*aij->ops->restorearray)(A,array)); 4427 } else { 4428 *array = NULL; 4429 } 4430 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4431 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4432 PetscFunctionReturn(0); 4433 } 4434 4435 /*@C 4436 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored 4437 4438 Not Collective 4439 4440 Input Parameter: 4441 . mat - a MATSEQAIJ matrix 4442 4443 Output Parameter: 4444 . array - pointer to the data 4445 4446 Level: intermediate 4447 4448 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4449 @*/ 4450 PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array) 4451 { 4452 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4453 4454 PetscFunctionBegin; 4455 if (aij->ops->getarrayread) { 4456 PetscCall((*aij->ops->getarrayread)(A,array)); 4457 } else { 4458 *array = aij->a; 4459 } 4460 PetscFunctionReturn(0); 4461 } 4462 4463 /*@C 4464 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4465 4466 Not Collective 4467 4468 Input Parameter: 4469 . mat - a MATSEQAIJ matrix 4470 4471 Output Parameter: 4472 . array - pointer to the data 4473 4474 Level: intermediate 4475 4476 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4477 @*/ 4478 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array) 4479 { 4480 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4481 4482 PetscFunctionBegin; 4483 if (aij->ops->restorearrayread) { 4484 PetscCall((*aij->ops->restorearrayread)(A,array)); 4485 } else { 4486 *array = NULL; 4487 } 4488 PetscFunctionReturn(0); 4489 } 4490 4491 /*@C 4492 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a MATSEQAIJ matrix is stored 4493 4494 Not Collective 4495 4496 Input Parameter: 4497 . mat - a MATSEQAIJ matrix 4498 4499 Output Parameter: 4500 . array - pointer to the data 4501 4502 Level: intermediate 4503 4504 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4505 @*/ 4506 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A,PetscScalar **array) 4507 { 4508 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4509 4510 PetscFunctionBegin; 4511 if (aij->ops->getarraywrite) { 4512 PetscCall((*aij->ops->getarraywrite)(A,array)); 4513 } else { 4514 *array = aij->a; 4515 } 4516 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4517 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4518 PetscFunctionReturn(0); 4519 } 4520 4521 /*@C 4522 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4523 4524 Not Collective 4525 4526 Input Parameter: 4527 . mat - a MATSEQAIJ matrix 4528 4529 Output Parameter: 4530 . array - pointer to the data 4531 4532 Level: intermediate 4533 4534 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4535 @*/ 4536 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A,PetscScalar **array) 4537 { 4538 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4539 4540 PetscFunctionBegin; 4541 if (aij->ops->restorearraywrite) { 4542 PetscCall((*aij->ops->restorearraywrite)(A,array)); 4543 } else { 4544 *array = NULL; 4545 } 4546 PetscFunctionReturn(0); 4547 } 4548 4549 /*@C 4550 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the SEQAIJ matrix 4551 4552 Not Collective 4553 4554 Input Parameter: 4555 . mat - a matrix of type MATSEQAIJ or its subclasses 4556 4557 Output Parameters: 4558 + i - row map array of the matrix 4559 . j - column index array of the matrix 4560 . a - data array of the matrix 4561 - memtype - memory type of the arrays 4562 4563 Notes: 4564 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4565 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4566 4567 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4568 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4569 4570 Level: Developer 4571 4572 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4573 @*/ 4574 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat,const PetscInt **i,const PetscInt **j,PetscScalar **a,PetscMemType *mtype) 4575 { 4576 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 4577 4578 PetscFunctionBegin; 4579 PetscCheck(mat->preallocated,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"matrix is not preallocated"); 4580 if (aij->ops->getcsrandmemtype) { 4581 PetscCall((*aij->ops->getcsrandmemtype)(mat,i,j,a,mtype)); 4582 } else { 4583 if (i) *i = aij->i; 4584 if (j) *j = aij->j; 4585 if (a) *a = aij->a; 4586 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4587 } 4588 PetscFunctionReturn(0); 4589 } 4590 4591 /*@C 4592 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4593 4594 Not Collective 4595 4596 Input Parameter: 4597 . mat - a MATSEQAIJ matrix 4598 4599 Output Parameter: 4600 . nz - the maximum number of nonzeros in any row 4601 4602 Level: intermediate 4603 4604 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4605 @*/ 4606 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4607 { 4608 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4609 4610 PetscFunctionBegin; 4611 *nz = aij->rmax; 4612 PetscFunctionReturn(0); 4613 } 4614 4615 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[]) 4616 { 4617 MPI_Comm comm; 4618 PetscInt *i,*j; 4619 PetscInt M,N,row; 4620 PetscCount k,p,q,nneg,nnz,start,end; /* Index the coo array, so use PetscCount as their type */ 4621 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4622 PetscInt *Aj; 4623 PetscScalar *Aa; 4624 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ*)(mat->data); 4625 MatType rtype; 4626 PetscCount *perm,*jmap; 4627 4628 PetscFunctionBegin; 4629 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4630 PetscCall(PetscObjectGetComm((PetscObject)mat,&comm)); 4631 PetscCall(MatGetSize(mat,&M,&N)); 4632 i = coo_i; 4633 j = coo_j; 4634 PetscCall(PetscMalloc1(coo_n,&perm)); 4635 for (k=0; k<coo_n; k++) { /* Ignore entries with negative row or col indices */ 4636 if (j[k] < 0) i[k] = -1; 4637 perm[k] = k; 4638 } 4639 4640 /* Sort by row */ 4641 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n,i,j,perm)); 4642 for (k=0; k<coo_n; k++) {if (i[k] >= 0) break;} /* Advance k to the first row with a non-negative index */ 4643 nneg = k; 4644 PetscCall(PetscMalloc1(coo_n-nneg+1,&jmap)); /* +1 to make a CSR-like data structure. jmap[i] originally is the number of repeats for i-th nonzero */ 4645 nnz = 0; /* Total number of unique nonzeros to be counted */ 4646 jmap++; /* Inc jmap by 1 for convinience */ 4647 4648 PetscCall(PetscCalloc1(M+1,&Ai)); /* CSR of A */ 4649 PetscCall(PetscMalloc1(coo_n-nneg,&Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4650 4651 /* In each row, sort by column, then unique column indices to get row length */ 4652 Ai++; /* Inc by 1 for convinience */ 4653 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4654 while (k<coo_n) { 4655 row = i[k]; 4656 start = k; /* [start,end) indices for this row */ 4657 while (k<coo_n && i[k] == row) k++; 4658 end = k; 4659 PetscCall(PetscSortIntWithCountArray(end-start,j+start,perm+start)); 4660 /* Find number of unique col entries in this row */ 4661 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4662 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4663 Ai[row] = 1; 4664 nnz++; 4665 4666 for (p=start+1; p<end; p++) { /* Scan remaining nonzero in this row */ 4667 if (j[p] != j[p-1]) { /* Meet a new nonzero */ 4668 q++; 4669 jmap[q] = 1; 4670 Aj[q] = j[p]; 4671 Ai[row]++; 4672 nnz++; 4673 } else { 4674 jmap[q]++; 4675 } 4676 } 4677 q++; /* Move to next row and thus next unique nonzero */ 4678 } 4679 4680 Ai--; /* Back to the beginning of Ai[] */ 4681 for (k=0; k<M; k++) Ai[k+1] += Ai[k]; 4682 jmap--; /* Back to the beginning of jmap[] */ 4683 jmap[0] = 0; 4684 for (k=0; k<nnz; k++) jmap[k+1] += jmap[k]; 4685 if (nnz < coo_n-nneg) { /* Realloc with actual number of unique nonzeros */ 4686 PetscCount *jmap_new; 4687 PetscInt *Aj_new; 4688 4689 PetscCall(PetscMalloc1(nnz+1,&jmap_new)); 4690 PetscCall(PetscArraycpy(jmap_new,jmap,nnz+1)); 4691 PetscCall(PetscFree(jmap)); 4692 jmap = jmap_new; 4693 4694 PetscCall(PetscMalloc1(nnz,&Aj_new)); 4695 PetscCall(PetscArraycpy(Aj_new,Aj,nnz)); 4696 PetscCall(PetscFree(Aj)); 4697 Aj = Aj_new; 4698 } 4699 4700 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4701 PetscCount *perm_new; 4702 4703 PetscCall(PetscMalloc1(coo_n-nneg,&perm_new)); 4704 PetscCall(PetscArraycpy(perm_new,perm+nneg,coo_n-nneg)); 4705 PetscCall(PetscFree(perm)); 4706 perm = perm_new; 4707 } 4708 4709 PetscCall(MatGetRootType_Private(mat,&rtype)); 4710 PetscCall(PetscCalloc1(nnz,&Aa)); /* Zero the matrix */ 4711 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF,M,N,Ai,Aj,Aa,rtype,mat)); 4712 4713 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4714 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4715 /* Record COO fields */ 4716 seqaij->coo_n = coo_n; 4717 seqaij->Atot = coo_n-nneg; /* Annz is seqaij->nz, so no need to record that again */ 4718 seqaij->jmap = jmap; /* of length nnz+1 */ 4719 seqaij->perm = perm; 4720 PetscFunctionReturn(0); 4721 } 4722 4723 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A,const PetscScalar v[],InsertMode imode) 4724 { 4725 Mat_SeqAIJ *aseq = (Mat_SeqAIJ*)A->data; 4726 PetscCount i,j,Annz = aseq->nz; 4727 PetscCount *perm = aseq->perm,*jmap = aseq->jmap; 4728 PetscScalar *Aa; 4729 4730 PetscFunctionBegin; 4731 PetscCall(MatSeqAIJGetArray(A,&Aa)); 4732 for (i=0; i<Annz; i++) { 4733 PetscScalar sum = 0.0; 4734 for (j=jmap[i]; j<jmap[i+1]; j++) sum += v[perm[j]]; 4735 Aa[i] = (imode == INSERT_VALUES? 0.0 : Aa[i]) + sum; 4736 } 4737 PetscCall(MatSeqAIJRestoreArray(A,&Aa)); 4738 PetscFunctionReturn(0); 4739 } 4740 4741 #if defined(PETSC_HAVE_CUDA) 4742 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*); 4743 #endif 4744 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4745 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*); 4746 #endif 4747 4748 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4749 { 4750 Mat_SeqAIJ *b; 4751 PetscMPIInt size; 4752 4753 PetscFunctionBegin; 4754 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 4755 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4756 4757 PetscCall(PetscNewLog(B,&b)); 4758 4759 B->data = (void*)b; 4760 4761 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 4762 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4763 4764 b->row = NULL; 4765 b->col = NULL; 4766 b->icol = NULL; 4767 b->reallocs = 0; 4768 b->ignorezeroentries = PETSC_FALSE; 4769 b->roworiented = PETSC_TRUE; 4770 b->nonew = 0; 4771 b->diag = NULL; 4772 b->solve_work = NULL; 4773 B->spptr = NULL; 4774 b->saved_values = NULL; 4775 b->idiag = NULL; 4776 b->mdiag = NULL; 4777 b->ssor_work = NULL; 4778 b->omega = 1.0; 4779 b->fshift = 0.0; 4780 b->idiagvalid = PETSC_FALSE; 4781 b->ibdiagvalid = PETSC_FALSE; 4782 b->keepnonzeropattern = PETSC_FALSE; 4783 4784 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4785 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4786 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ)); 4787 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ)); 4788 #endif 4789 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ)); 4790 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ)); 4791 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ)); 4792 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ)); 4793 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ)); 4794 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM)); 4795 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL)); 4796 #if defined(PETSC_HAVE_MKL_SPARSE) 4797 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL)); 4798 #endif 4799 #if defined(PETSC_HAVE_CUDA) 4800 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4801 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4802 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ)); 4803 #endif 4804 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4805 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos)); 4806 #endif 4807 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL)); 4808 #if defined(PETSC_HAVE_ELEMENTAL) 4809 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental)); 4810 #endif 4811 #if defined(PETSC_HAVE_SCALAPACK) 4812 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK)); 4813 #endif 4814 #if defined(PETSC_HAVE_HYPRE) 4815 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE)); 4816 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4817 #endif 4818 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense)); 4819 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL)); 4820 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS)); 4821 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ)); 4822 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ)); 4823 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ)); 4824 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ)); 4825 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4826 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ)); 4827 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ)); 4828 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ)); 4829 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4830 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ)); 4831 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_SeqAIJ)); 4832 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_SeqAIJ)); 4833 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4834 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4835 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4836 PetscFunctionReturn(0); 4837 } 4838 4839 /* 4840 Given a matrix generated with MatGetFactor() duplicates all the information in A into C 4841 */ 4842 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4843 { 4844 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data; 4845 PetscInt m = A->rmap->n,i; 4846 4847 PetscFunctionBegin; 4848 PetscCheck(A->assembled || cpvalues == MAT_DO_NOT_COPY_VALUES,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot duplicate unassembled matrix"); 4849 4850 C->factortype = A->factortype; 4851 c->row = NULL; 4852 c->col = NULL; 4853 c->icol = NULL; 4854 c->reallocs = 0; 4855 4856 C->assembled = A->assembled; 4857 C->preallocated = A->preallocated; 4858 4859 if (A->preallocated) { 4860 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 4861 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 4862 4863 PetscCall(PetscMalloc1(m,&c->imax)); 4864 PetscCall(PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt))); 4865 PetscCall(PetscMalloc1(m,&c->ilen)); 4866 PetscCall(PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt))); 4867 PetscCall(PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt))); 4868 4869 /* allocate the matrix space */ 4870 if (mallocmatspace) { 4871 PetscCall(PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i)); 4872 PetscCall(PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt))); 4873 4874 c->singlemalloc = PETSC_TRUE; 4875 4876 PetscCall(PetscArraycpy(c->i,a->i,m+1)); 4877 if (m > 0) { 4878 PetscCall(PetscArraycpy(c->j,a->j,a->i[m])); 4879 if (cpvalues == MAT_COPY_VALUES) { 4880 const PetscScalar *aa; 4881 4882 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4883 PetscCall(PetscArraycpy(c->a,aa,a->i[m])); 4884 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4885 } else { 4886 PetscCall(PetscArrayzero(c->a,a->i[m])); 4887 } 4888 } 4889 } 4890 4891 c->ignorezeroentries = a->ignorezeroentries; 4892 c->roworiented = a->roworiented; 4893 c->nonew = a->nonew; 4894 if (a->diag) { 4895 PetscCall(PetscMalloc1(m+1,&c->diag)); 4896 PetscCall(PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt))); 4897 PetscCall(PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt))); 4898 } else c->diag = NULL; 4899 4900 c->solve_work = NULL; 4901 c->saved_values = NULL; 4902 c->idiag = NULL; 4903 c->ssor_work = NULL; 4904 c->keepnonzeropattern = a->keepnonzeropattern; 4905 c->free_a = PETSC_TRUE; 4906 c->free_ij = PETSC_TRUE; 4907 4908 c->rmax = a->rmax; 4909 c->nz = a->nz; 4910 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4911 4912 c->compressedrow.use = a->compressedrow.use; 4913 c->compressedrow.nrows = a->compressedrow.nrows; 4914 if (a->compressedrow.use) { 4915 i = a->compressedrow.nrows; 4916 PetscCall(PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex)); 4917 PetscCall(PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1)); 4918 PetscCall(PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i)); 4919 } else { 4920 c->compressedrow.use = PETSC_FALSE; 4921 c->compressedrow.i = NULL; 4922 c->compressedrow.rindex = NULL; 4923 } 4924 c->nonzerorowcnt = a->nonzerorowcnt; 4925 C->nonzerostate = A->nonzerostate; 4926 4927 PetscCall(MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C)); 4928 } 4929 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 4930 PetscFunctionReturn(0); 4931 } 4932 4933 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4934 { 4935 PetscFunctionBegin; 4936 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 4937 PetscCall(MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 4938 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4939 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 4940 } 4941 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 4942 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE)); 4943 PetscFunctionReturn(0); 4944 } 4945 4946 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4947 { 4948 PetscBool isbinary, ishdf5; 4949 4950 PetscFunctionBegin; 4951 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4952 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4953 /* force binary viewer to load .info file if it has not yet done so */ 4954 PetscCall(PetscViewerSetUp(viewer)); 4955 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 4956 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5)); 4957 if (isbinary) { 4958 PetscCall(MatLoad_SeqAIJ_Binary(newMat,viewer)); 4959 } else if (ishdf5) { 4960 #if defined(PETSC_HAVE_HDF5) 4961 PetscCall(MatLoad_AIJ_HDF5(newMat,viewer)); 4962 #else 4963 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4964 #endif 4965 } else { 4966 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 4967 } 4968 PetscFunctionReturn(0); 4969 } 4970 4971 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 4972 { 4973 Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data; 4974 PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i; 4975 4976 PetscFunctionBegin; 4977 PetscCall(PetscViewerSetUp(viewer)); 4978 4979 /* read in matrix header */ 4980 PetscCall(PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT)); 4981 PetscCheck(header[0] == MAT_FILE_CLASSID,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 4982 M = header[1]; N = header[2]; nz = header[3]; 4983 PetscCheck(M >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%" PetscInt_FMT ") in file is negative",M); 4984 PetscCheck(N >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%" PetscInt_FMT ") in file is negative",N); 4985 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 4986 4987 /* set block sizes from the viewer's .info file */ 4988 PetscCall(MatLoad_Binary_BlockSizes(mat,viewer)); 4989 /* set local and global sizes if not set already */ 4990 if (mat->rmap->n < 0) mat->rmap->n = M; 4991 if (mat->cmap->n < 0) mat->cmap->n = N; 4992 if (mat->rmap->N < 0) mat->rmap->N = M; 4993 if (mat->cmap->N < 0) mat->cmap->N = N; 4994 PetscCall(PetscLayoutSetUp(mat->rmap)); 4995 PetscCall(PetscLayoutSetUp(mat->cmap)); 4996 4997 /* check if the matrix sizes are correct */ 4998 PetscCall(MatGetSize(mat,&rows,&cols)); 4999 PetscCheck(M == rows && N == cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5000 5001 /* read in row lengths */ 5002 PetscCall(PetscMalloc1(M,&rowlens)); 5003 PetscCall(PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT)); 5004 /* check if sum(rowlens) is same as nz */ 5005 sum = 0; for (i=0; i<M; i++) sum += rowlens[i]; 5006 PetscCheck(sum == nz,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT,nz,sum); 5007 /* preallocate and check sizes */ 5008 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens)); 5009 PetscCall(MatGetSize(mat,&rows,&cols)); 5010 PetscCheck(M == rows && N == cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5011 /* store row lengths */ 5012 PetscCall(PetscArraycpy(a->ilen,rowlens,M)); 5013 PetscCall(PetscFree(rowlens)); 5014 5015 /* fill in "i" row pointers */ 5016 a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i]; 5017 /* read in "j" column indices */ 5018 PetscCall(PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT)); 5019 /* read in "a" nonzero values */ 5020 PetscCall(PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR)); 5021 5022 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 5023 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 5024 PetscFunctionReturn(0); 5025 } 5026 5027 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 5028 { 5029 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 5030 const PetscScalar *aa,*ba; 5031 #if defined(PETSC_USE_COMPLEX) 5032 PetscInt k; 5033 #endif 5034 5035 PetscFunctionBegin; 5036 /* If the matrix dimensions are not equal,or no of nonzeros */ 5037 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 5038 *flg = PETSC_FALSE; 5039 PetscFunctionReturn(0); 5040 } 5041 5042 /* if the a->i are the same */ 5043 PetscCall(PetscArraycmp(a->i,b->i,A->rmap->n+1,flg)); 5044 if (!*flg) PetscFunctionReturn(0); 5045 5046 /* if a->j are the same */ 5047 PetscCall(PetscArraycmp(a->j,b->j,a->nz,flg)); 5048 if (!*flg) PetscFunctionReturn(0); 5049 5050 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 5051 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 5052 /* if a->a are the same */ 5053 #if defined(PETSC_USE_COMPLEX) 5054 for (k=0; k<a->nz; k++) { 5055 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5056 *flg = PETSC_FALSE; 5057 PetscFunctionReturn(0); 5058 } 5059 } 5060 #else 5061 PetscCall(PetscArraycmp(aa,ba,a->nz,flg)); 5062 #endif 5063 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 5064 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 5065 PetscFunctionReturn(0); 5066 } 5067 5068 /*@ 5069 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 5070 provided by the user. 5071 5072 Collective 5073 5074 Input Parameters: 5075 + comm - must be an MPI communicator of size 1 5076 . m - number of rows 5077 . n - number of columns 5078 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5079 . j - column indices 5080 - a - matrix values 5081 5082 Output Parameter: 5083 . mat - the matrix 5084 5085 Level: intermediate 5086 5087 Notes: 5088 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5089 once the matrix is destroyed and not before 5090 5091 You cannot set new nonzero locations into this matrix, that will generate an error. 5092 5093 The i and j indices are 0 based 5094 5095 The format which is used for the sparse matrix input, is equivalent to a 5096 row-major ordering.. i.e for the following matrix, the input data expected is 5097 as shown 5098 5099 $ 1 0 0 5100 $ 2 0 3 5101 $ 4 5 6 5102 $ 5103 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5104 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5105 $ v = {1,2,3,4,5,6} [size = 6] 5106 5107 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateMPIAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()` 5108 5109 @*/ 5110 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 5111 { 5112 PetscInt ii; 5113 Mat_SeqAIJ *aij; 5114 PetscInt jj; 5115 5116 PetscFunctionBegin; 5117 PetscCheck(m <= 0 || i[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 5118 PetscCall(MatCreate(comm,mat)); 5119 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5120 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5121 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5122 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL)); 5123 aij = (Mat_SeqAIJ*)(*mat)->data; 5124 PetscCall(PetscMalloc1(m,&aij->imax)); 5125 PetscCall(PetscMalloc1(m,&aij->ilen)); 5126 5127 aij->i = i; 5128 aij->j = j; 5129 aij->a = a; 5130 aij->singlemalloc = PETSC_FALSE; 5131 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5132 aij->free_a = PETSC_FALSE; 5133 aij->free_ij = PETSC_FALSE; 5134 5135 for (ii=0,aij->nonzerorowcnt=0,aij->rmax=0; ii<m; ii++) { 5136 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 5137 if (PetscDefined(USE_DEBUG)) { 5138 PetscCheck(i[ii+1] - i[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT,ii,i[ii+1] - i[ii]); 5139 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 5140 PetscCheck(j[jj] >= j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is not sorted",jj-i[ii],j[jj],ii); 5141 PetscCheck(j[jj] != j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is identical to previous entry",jj-i[ii],j[jj],ii); 5142 } 5143 } 5144 } 5145 if (PetscDefined(USE_DEBUG)) { 5146 for (ii=0; ii<aij->i[m]; ii++) { 5147 PetscCheck(j[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5148 PetscCheck(j[ii] <= n - 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5149 } 5150 } 5151 5152 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5153 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5154 PetscFunctionReturn(0); 5155 } 5156 5157 /*@ 5158 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 5159 provided by the user. 5160 5161 Collective 5162 5163 Input Parameters: 5164 + comm - must be an MPI communicator of size 1 5165 . m - number of rows 5166 . n - number of columns 5167 . i - row indices 5168 . j - column indices 5169 . a - matrix values 5170 . nz - number of nonzeros 5171 - idx - if the i and j indices start with 1 use PETSC_TRUE otherwise use PETSC_FALSE 5172 5173 Output Parameter: 5174 . mat - the matrix 5175 5176 Level: intermediate 5177 5178 Example: 5179 For the following matrix, the input data expected is as shown (using 0 based indexing) 5180 .vb 5181 1 0 0 5182 2 0 3 5183 4 5 6 5184 5185 i = {0,1,1,2,2,2} 5186 j = {0,0,2,0,1,2} 5187 v = {1,2,3,4,5,6} 5188 .ve 5189 5190 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateSeqAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()`, `MatSetValuesCOO()` 5191 5192 @*/ 5193 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 5194 { 5195 PetscInt ii, *nnz, one = 1,row,col; 5196 5197 PetscFunctionBegin; 5198 PetscCall(PetscCalloc1(m,&nnz)); 5199 for (ii = 0; ii < nz; ii++) { 5200 nnz[i[ii] - !!idx] += 1; 5201 } 5202 PetscCall(MatCreate(comm,mat)); 5203 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5204 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5205 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz)); 5206 for (ii = 0; ii < nz; ii++) { 5207 if (idx) { 5208 row = i[ii] - 1; 5209 col = j[ii] - 1; 5210 } else { 5211 row = i[ii]; 5212 col = j[ii]; 5213 } 5214 PetscCall(MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES)); 5215 } 5216 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5217 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5218 PetscCall(PetscFree(nnz)); 5219 PetscFunctionReturn(0); 5220 } 5221 5222 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5223 { 5224 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 5225 5226 PetscFunctionBegin; 5227 a->idiagvalid = PETSC_FALSE; 5228 a->ibdiagvalid = PETSC_FALSE; 5229 5230 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5231 PetscFunctionReturn(0); 5232 } 5233 5234 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 5235 { 5236 PetscFunctionBegin; 5237 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat)); 5238 PetscFunctionReturn(0); 5239 } 5240 5241 /* 5242 Permute A into C's *local* index space using rowemb,colemb. 5243 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5244 of [0,m), colemb is in [0,n). 5245 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5246 */ 5247 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 5248 { 5249 /* If making this function public, change the error returned in this function away from _PLIB. */ 5250 Mat_SeqAIJ *Baij; 5251 PetscBool seqaij; 5252 PetscInt m,n,*nz,i,j,count; 5253 PetscScalar v; 5254 const PetscInt *rowindices,*colindices; 5255 5256 PetscFunctionBegin; 5257 if (!B) PetscFunctionReturn(0); 5258 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5259 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij)); 5260 PetscCheck(seqaij,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 5261 if (rowemb) { 5262 PetscCall(ISGetLocalSize(rowemb,&m)); 5263 PetscCheck(m == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row IS of size %" PetscInt_FMT " is incompatible with matrix row size %" PetscInt_FMT,m,B->rmap->n); 5264 } else { 5265 PetscCheck(C->rmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 5266 } 5267 if (colemb) { 5268 PetscCall(ISGetLocalSize(colemb,&n)); 5269 PetscCheck(n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diag col IS of size %" PetscInt_FMT " is incompatible with input matrix col size %" PetscInt_FMT,n,B->cmap->n); 5270 } else { 5271 PetscCheck(C->cmap->n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 5272 } 5273 5274 Baij = (Mat_SeqAIJ*)(B->data); 5275 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5276 PetscCall(PetscMalloc1(B->rmap->n,&nz)); 5277 for (i=0; i<B->rmap->n; i++) { 5278 nz[i] = Baij->i[i+1] - Baij->i[i]; 5279 } 5280 PetscCall(MatSeqAIJSetPreallocation(C,0,nz)); 5281 PetscCall(PetscFree(nz)); 5282 } 5283 if (pattern == SUBSET_NONZERO_PATTERN) { 5284 PetscCall(MatZeroEntries(C)); 5285 } 5286 count = 0; 5287 rowindices = NULL; 5288 colindices = NULL; 5289 if (rowemb) { 5290 PetscCall(ISGetIndices(rowemb,&rowindices)); 5291 } 5292 if (colemb) { 5293 PetscCall(ISGetIndices(colemb,&colindices)); 5294 } 5295 for (i=0; i<B->rmap->n; i++) { 5296 PetscInt row; 5297 row = i; 5298 if (rowindices) row = rowindices[i]; 5299 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 5300 PetscInt col; 5301 col = Baij->j[count]; 5302 if (colindices) col = colindices[col]; 5303 v = Baij->a[count]; 5304 PetscCall(MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES)); 5305 ++count; 5306 } 5307 } 5308 /* FIXME: set C's nonzerostate correctly. */ 5309 /* Assembly for C is necessary. */ 5310 C->preallocated = PETSC_TRUE; 5311 C->assembled = PETSC_TRUE; 5312 C->was_assembled = PETSC_FALSE; 5313 PetscFunctionReturn(0); 5314 } 5315 5316 PetscFunctionList MatSeqAIJList = NULL; 5317 5318 /*@C 5319 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 5320 5321 Collective on Mat 5322 5323 Input Parameters: 5324 + mat - the matrix object 5325 - matype - matrix type 5326 5327 Options Database Key: 5328 . -mat_seqai_type <method> - for example seqaijcrl 5329 5330 Level: intermediate 5331 5332 .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat` 5333 @*/ 5334 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5335 { 5336 PetscBool sametype; 5337 PetscErrorCode (*r)(Mat,MatType,MatReuse,Mat*); 5338 5339 PetscFunctionBegin; 5340 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 5341 PetscCall(PetscObjectTypeCompare((PetscObject)mat,matype,&sametype)); 5342 if (sametype) PetscFunctionReturn(0); 5343 5344 PetscCall(PetscFunctionListFind(MatSeqAIJList,matype,&r)); 5345 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 5346 PetscCall((*r)(mat,matype,MAT_INPLACE_MATRIX,&mat)); 5347 PetscFunctionReturn(0); 5348 } 5349 5350 /*@C 5351 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 5352 5353 Not Collective 5354 5355 Input Parameters: 5356 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 5357 - function - routine to convert to subtype 5358 5359 Notes: 5360 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 5361 5362 Then, your matrix can be chosen with the procedural interface at runtime via the option 5363 $ -mat_seqaij_type my_mat 5364 5365 Level: advanced 5366 5367 .seealso: `MatSeqAIJRegisterAll()` 5368 5369 Level: advanced 5370 @*/ 5371 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 5372 { 5373 PetscFunctionBegin; 5374 PetscCall(MatInitializePackage()); 5375 PetscCall(PetscFunctionListAdd(&MatSeqAIJList,sname,function)); 5376 PetscFunctionReturn(0); 5377 } 5378 5379 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5380 5381 /*@C 5382 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 5383 5384 Not Collective 5385 5386 Level: advanced 5387 5388 .seealso: `MatRegisterAll()`, `MatSeqAIJRegister()` 5389 @*/ 5390 PetscErrorCode MatSeqAIJRegisterAll(void) 5391 { 5392 PetscFunctionBegin; 5393 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5394 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5395 5396 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5397 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5398 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5399 #if defined(PETSC_HAVE_MKL_SPARSE) 5400 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5401 #endif 5402 #if defined(PETSC_HAVE_CUDA) 5403 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5404 #endif 5405 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5406 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5407 #endif 5408 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5409 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5410 #endif 5411 PetscFunctionReturn(0); 5412 } 5413 5414 /* 5415 Special version for direct calls from Fortran 5416 */ 5417 #include <petsc/private/fortranimpl.h> 5418 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5419 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5420 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5421 #define matsetvaluesseqaij_ matsetvaluesseqaij 5422 #endif 5423 5424 /* Change these macros so can be used in void function */ 5425 5426 /* Change these macros so can be used in void function */ 5427 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5428 #undef PetscCall 5429 #define PetscCall(...) do { \ 5430 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5431 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5432 *_ierr = PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr_msv_mpiaij,PETSC_ERROR_REPEAT," "); \ 5433 return; \ 5434 } \ 5435 } while (0) 5436 5437 #undef SETERRQ 5438 #define SETERRQ(comm,ierr,...) do { \ 5439 *_ierr = PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__); \ 5440 return; \ 5441 } while (0) 5442 5443 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 5444 { 5445 Mat A = *AA; 5446 PetscInt m = *mm, n = *nn; 5447 InsertMode is = *isis; 5448 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 5449 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 5450 PetscInt *imax,*ai,*ailen; 5451 PetscInt *aj,nonew = a->nonew,lastcol = -1; 5452 MatScalar *ap,value,*aa; 5453 PetscBool ignorezeroentries = a->ignorezeroentries; 5454 PetscBool roworiented = a->roworiented; 5455 5456 PetscFunctionBegin; 5457 MatCheckPreallocated(A,1); 5458 imax = a->imax; 5459 ai = a->i; 5460 ailen = a->ilen; 5461 aj = a->j; 5462 aa = a->a; 5463 5464 for (k=0; k<m; k++) { /* loop over added rows */ 5465 row = im[k]; 5466 if (row < 0) continue; 5467 PetscCheck(row < A->rmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 5468 rp = aj + ai[row]; ap = aa + ai[row]; 5469 rmax = imax[row]; nrow = ailen[row]; 5470 low = 0; 5471 high = nrow; 5472 for (l=0; l<n; l++) { /* loop over added columns */ 5473 if (in[l] < 0) continue; 5474 PetscCheck(in[l] < A->cmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 5475 col = in[l]; 5476 if (roworiented) value = v[l + k*n]; 5477 else value = v[k + l*m]; 5478 5479 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5480 5481 if (col <= lastcol) low = 0; 5482 else high = nrow; 5483 lastcol = col; 5484 while (high-low > 5) { 5485 t = (low+high)/2; 5486 if (rp[t] > col) high = t; 5487 else low = t; 5488 } 5489 for (i=low; i<high; i++) { 5490 if (rp[i] > col) break; 5491 if (rp[i] == col) { 5492 if (is == ADD_VALUES) ap[i] += value; 5493 else ap[i] = value; 5494 goto noinsert; 5495 } 5496 } 5497 if (value == 0.0 && ignorezeroentries) goto noinsert; 5498 if (nonew == 1) goto noinsert; 5499 PetscCheck(nonew != -1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 5500 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 5501 N = nrow++ - 1; a->nz++; high++; 5502 /* shift up all the later entries in this row */ 5503 for (ii=N; ii>=i; ii--) { 5504 rp[ii+1] = rp[ii]; 5505 ap[ii+1] = ap[ii]; 5506 } 5507 rp[i] = col; 5508 ap[i] = value; 5509 A->nonzerostate++; 5510 noinsert:; 5511 low = i + 1; 5512 } 5513 ailen[row] = nrow; 5514 } 5515 PetscFunctionReturnVoid(); 5516 } 5517 /* Undefining these here since they were redefined from their original definition above! No 5518 * other PETSc functions should be defined past this point, as it is impossible to recover the 5519 * original definitions */ 5520 #undef PetscCall 5521 #undef SETERRQ 5522