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