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