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