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