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