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,"MatIsTranspose_C",NULL);CHKERRQ(ierr); 1125 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 1126 ierr = PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL);CHKERRQ(ierr); 1127 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 1128 ierr = PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);CHKERRQ(ierr); 1129 PetscFunctionReturn(0); 1130 } 1131 1132 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1133 { 1134 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1135 PetscErrorCode ierr; 1136 1137 PetscFunctionBegin; 1138 switch (op) { 1139 case MAT_ROW_ORIENTED: 1140 a->roworiented = flg; 1141 break; 1142 case MAT_KEEP_NONZERO_PATTERN: 1143 a->keepnonzeropattern = flg; 1144 break; 1145 case MAT_NEW_NONZERO_LOCATIONS: 1146 a->nonew = (flg ? 0 : 1); 1147 break; 1148 case MAT_NEW_NONZERO_LOCATION_ERR: 1149 a->nonew = (flg ? -1 : 0); 1150 break; 1151 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1152 a->nonew = (flg ? -2 : 0); 1153 break; 1154 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1155 a->nounused = (flg ? -1 : 0); 1156 break; 1157 case MAT_IGNORE_ZERO_ENTRIES: 1158 a->ignorezeroentries = flg; 1159 break; 1160 case MAT_SPD: 1161 case MAT_SYMMETRIC: 1162 case MAT_STRUCTURALLY_SYMMETRIC: 1163 case MAT_HERMITIAN: 1164 case MAT_SYMMETRY_ETERNAL: 1165 case MAT_STRUCTURE_ONLY: 1166 /* These options are handled directly by MatSetOption() */ 1167 break; 1168 case MAT_NEW_DIAGONALS: 1169 case MAT_IGNORE_OFF_PROC_ENTRIES: 1170 case MAT_USE_HASH_TABLE: 1171 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1172 break; 1173 case MAT_USE_INODES: 1174 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 1175 break; 1176 case MAT_SUBMAT_SINGLEIS: 1177 A->submat_singleis = flg; 1178 break; 1179 default: 1180 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1181 } 1182 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 1183 PetscFunctionReturn(0); 1184 } 1185 1186 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1187 { 1188 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1189 PetscErrorCode ierr; 1190 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 1191 PetscScalar *aa=a->a,*x,zero=0.0; 1192 1193 PetscFunctionBegin; 1194 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1195 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1196 1197 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1198 PetscInt *diag=a->diag; 1199 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1200 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1201 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1202 PetscFunctionReturn(0); 1203 } 1204 1205 ierr = VecSet(v,zero);CHKERRQ(ierr); 1206 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1207 for (i=0; i<n; i++) { 1208 nz = ai[i+1] - ai[i]; 1209 if (!nz) x[i] = 0.0; 1210 for (j=ai[i]; j<ai[i+1]; j++) { 1211 if (aj[j] == i) { 1212 x[i] = aa[j]; 1213 break; 1214 } 1215 } 1216 } 1217 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1218 PetscFunctionReturn(0); 1219 } 1220 1221 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1222 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1223 { 1224 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1225 PetscScalar *y; 1226 const PetscScalar *x; 1227 PetscErrorCode ierr; 1228 PetscInt m = A->rmap->n; 1229 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1230 const MatScalar *v; 1231 PetscScalar alpha; 1232 PetscInt n,i,j; 1233 const PetscInt *idx,*ii,*ridx=NULL; 1234 Mat_CompressedRow cprow = a->compressedrow; 1235 PetscBool usecprow = cprow.use; 1236 #endif 1237 1238 PetscFunctionBegin; 1239 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 1240 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1241 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1242 1243 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1244 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 1245 #else 1246 if (usecprow) { 1247 m = cprow.nrows; 1248 ii = cprow.i; 1249 ridx = cprow.rindex; 1250 } else { 1251 ii = a->i; 1252 } 1253 for (i=0; i<m; i++) { 1254 idx = a->j + ii[i]; 1255 v = a->a + ii[i]; 1256 n = ii[i+1] - ii[i]; 1257 if (usecprow) { 1258 alpha = x[ridx[i]]; 1259 } else { 1260 alpha = x[i]; 1261 } 1262 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1263 } 1264 #endif 1265 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1266 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1267 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1268 PetscFunctionReturn(0); 1269 } 1270 1271 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1272 { 1273 PetscErrorCode ierr; 1274 1275 PetscFunctionBegin; 1276 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 1277 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 1278 PetscFunctionReturn(0); 1279 } 1280 1281 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1282 1283 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1284 { 1285 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1286 PetscScalar *y; 1287 const PetscScalar *x; 1288 const MatScalar *aa; 1289 PetscErrorCode ierr; 1290 PetscInt m=A->rmap->n; 1291 const PetscInt *aj,*ii,*ridx=NULL; 1292 PetscInt n,i; 1293 PetscScalar sum; 1294 PetscBool usecprow=a->compressedrow.use; 1295 1296 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1297 #pragma disjoint(*x,*y,*aa) 1298 #endif 1299 1300 PetscFunctionBegin; 1301 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1302 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1303 ii = a->i; 1304 if (usecprow) { /* use compressed row format */ 1305 ierr = PetscMemzero(y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1306 m = a->compressedrow.nrows; 1307 ii = a->compressedrow.i; 1308 ridx = a->compressedrow.rindex; 1309 for (i=0; i<m; i++) { 1310 n = ii[i+1] - ii[i]; 1311 aj = a->j + ii[i]; 1312 aa = a->a + ii[i]; 1313 sum = 0.0; 1314 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1315 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1316 y[*ridx++] = sum; 1317 } 1318 } else { /* do not use compressed row format */ 1319 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1320 aj = a->j; 1321 aa = a->a; 1322 fortranmultaij_(&m,x,ii,aj,aa,y); 1323 #else 1324 for (i=0; i<m; i++) { 1325 n = ii[i+1] - ii[i]; 1326 aj = a->j + ii[i]; 1327 aa = a->a + ii[i]; 1328 sum = 0.0; 1329 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1330 y[i] = sum; 1331 } 1332 #endif 1333 } 1334 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1335 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1336 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1337 PetscFunctionReturn(0); 1338 } 1339 1340 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1341 { 1342 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1343 PetscScalar *y; 1344 const PetscScalar *x; 1345 const MatScalar *aa; 1346 PetscErrorCode ierr; 1347 PetscInt m=A->rmap->n; 1348 const PetscInt *aj,*ii,*ridx=NULL; 1349 PetscInt n,i,nonzerorow=0; 1350 PetscScalar sum; 1351 PetscBool usecprow=a->compressedrow.use; 1352 1353 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1354 #pragma disjoint(*x,*y,*aa) 1355 #endif 1356 1357 PetscFunctionBegin; 1358 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1359 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1360 if (usecprow) { /* use compressed row format */ 1361 m = a->compressedrow.nrows; 1362 ii = a->compressedrow.i; 1363 ridx = a->compressedrow.rindex; 1364 for (i=0; i<m; i++) { 1365 n = ii[i+1] - ii[i]; 1366 aj = a->j + ii[i]; 1367 aa = a->a + ii[i]; 1368 sum = 0.0; 1369 nonzerorow += (n>0); 1370 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1371 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1372 y[*ridx++] = sum; 1373 } 1374 } else { /* do not use compressed row format */ 1375 ii = a->i; 1376 for (i=0; i<m; i++) { 1377 n = ii[i+1] - ii[i]; 1378 aj = a->j + ii[i]; 1379 aa = a->a + ii[i]; 1380 sum = 0.0; 1381 nonzerorow += (n>0); 1382 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1383 y[i] = sum; 1384 } 1385 } 1386 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1387 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1388 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1389 PetscFunctionReturn(0); 1390 } 1391 1392 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1393 { 1394 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1395 PetscScalar *y,*z; 1396 const PetscScalar *x; 1397 const MatScalar *aa; 1398 PetscErrorCode ierr; 1399 PetscInt m = A->rmap->n,*aj,*ii; 1400 PetscInt n,i,*ridx=NULL; 1401 PetscScalar sum; 1402 PetscBool usecprow=a->compressedrow.use; 1403 1404 PetscFunctionBegin; 1405 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1406 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1407 if (usecprow) { /* use compressed row format */ 1408 if (zz != yy) { 1409 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1410 } 1411 m = a->compressedrow.nrows; 1412 ii = a->compressedrow.i; 1413 ridx = a->compressedrow.rindex; 1414 for (i=0; i<m; i++) { 1415 n = ii[i+1] - ii[i]; 1416 aj = a->j + ii[i]; 1417 aa = a->a + ii[i]; 1418 sum = y[*ridx]; 1419 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1420 z[*ridx++] = sum; 1421 } 1422 } else { /* do not use compressed row format */ 1423 ii = a->i; 1424 for (i=0; i<m; i++) { 1425 n = ii[i+1] - ii[i]; 1426 aj = a->j + ii[i]; 1427 aa = a->a + ii[i]; 1428 sum = y[i]; 1429 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1430 z[i] = sum; 1431 } 1432 } 1433 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1434 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1435 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438 1439 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1440 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1441 { 1442 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1443 PetscScalar *y,*z; 1444 const PetscScalar *x; 1445 const MatScalar *aa; 1446 PetscErrorCode ierr; 1447 const PetscInt *aj,*ii,*ridx=NULL; 1448 PetscInt m = A->rmap->n,n,i; 1449 PetscScalar sum; 1450 PetscBool usecprow=a->compressedrow.use; 1451 1452 PetscFunctionBegin; 1453 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1454 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1455 if (usecprow) { /* use compressed row format */ 1456 if (zz != yy) { 1457 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1458 } 1459 m = a->compressedrow.nrows; 1460 ii = a->compressedrow.i; 1461 ridx = a->compressedrow.rindex; 1462 for (i=0; i<m; i++) { 1463 n = ii[i+1] - ii[i]; 1464 aj = a->j + ii[i]; 1465 aa = a->a + ii[i]; 1466 sum = y[*ridx]; 1467 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1468 z[*ridx++] = sum; 1469 } 1470 } else { /* do not use compressed row format */ 1471 ii = a->i; 1472 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1473 aj = a->j; 1474 aa = a->a; 1475 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1476 #else 1477 for (i=0; i<m; i++) { 1478 n = ii[i+1] - ii[i]; 1479 aj = a->j + ii[i]; 1480 aa = a->a + ii[i]; 1481 sum = y[i]; 1482 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1483 z[i] = sum; 1484 } 1485 #endif 1486 } 1487 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1488 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1489 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1490 #if defined(PETSC_HAVE_CUSP) 1491 /* 1492 ierr = VecView(xx,0);CHKERRQ(ierr); 1493 ierr = VecView(zz,0);CHKERRQ(ierr); 1494 ierr = MatView(A,0);CHKERRQ(ierr); 1495 */ 1496 #endif 1497 PetscFunctionReturn(0); 1498 } 1499 1500 /* 1501 Adds diagonal pointers to sparse matrix structure. 1502 */ 1503 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1504 { 1505 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1506 PetscErrorCode ierr; 1507 PetscInt i,j,m = A->rmap->n; 1508 1509 PetscFunctionBegin; 1510 if (!a->diag) { 1511 ierr = PetscMalloc1(m,&a->diag);CHKERRQ(ierr); 1512 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1513 } 1514 for (i=0; i<A->rmap->n; i++) { 1515 a->diag[i] = a->i[i+1]; 1516 for (j=a->i[i]; j<a->i[i+1]; j++) { 1517 if (a->j[j] == i) { 1518 a->diag[i] = j; 1519 break; 1520 } 1521 } 1522 } 1523 PetscFunctionReturn(0); 1524 } 1525 1526 /* 1527 Checks for missing diagonals 1528 */ 1529 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1530 { 1531 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1532 PetscInt *diag,*ii = a->i,i; 1533 1534 PetscFunctionBegin; 1535 *missing = PETSC_FALSE; 1536 if (A->rmap->n > 0 && !ii) { 1537 *missing = PETSC_TRUE; 1538 if (d) *d = 0; 1539 PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n"); 1540 } else { 1541 diag = a->diag; 1542 for (i=0; i<A->rmap->n; i++) { 1543 if (diag[i] >= ii[i+1]) { 1544 *missing = PETSC_TRUE; 1545 if (d) *d = i; 1546 PetscInfo1(A,"Matrix is missing diagonal number %D\n",i); 1547 break; 1548 } 1549 } 1550 } 1551 PetscFunctionReturn(0); 1552 } 1553 1554 /* 1555 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1556 */ 1557 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1558 { 1559 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1560 PetscErrorCode ierr; 1561 PetscInt i,*diag,m = A->rmap->n; 1562 MatScalar *v = a->a; 1563 PetscScalar *idiag,*mdiag; 1564 1565 PetscFunctionBegin; 1566 if (a->idiagvalid) PetscFunctionReturn(0); 1567 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1568 diag = a->diag; 1569 if (!a->idiag) { 1570 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1571 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1572 v = a->a; 1573 } 1574 mdiag = a->mdiag; 1575 idiag = a->idiag; 1576 1577 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1578 for (i=0; i<m; i++) { 1579 mdiag[i] = v[diag[i]]; 1580 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1581 if (PetscRealPart(fshift)) { 1582 ierr = PetscInfo1(A,"Zero diagonal on row %D\n",i);CHKERRQ(ierr); 1583 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1584 A->factorerror_zeropivot_value = 0.0; 1585 A->factorerror_zeropivot_row = i; 1586 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1587 } 1588 idiag[i] = 1.0/v[diag[i]]; 1589 } 1590 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1591 } else { 1592 for (i=0; i<m; i++) { 1593 mdiag[i] = v[diag[i]]; 1594 idiag[i] = omega/(fshift + v[diag[i]]); 1595 } 1596 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1597 } 1598 a->idiagvalid = PETSC_TRUE; 1599 PetscFunctionReturn(0); 1600 } 1601 1602 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1603 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1604 { 1605 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1606 PetscScalar *x,d,sum,*t,scale; 1607 const MatScalar *v,*idiag=0,*mdiag; 1608 const PetscScalar *b, *bs,*xb, *ts; 1609 PetscErrorCode ierr; 1610 PetscInt n,m = A->rmap->n,i; 1611 const PetscInt *idx,*diag; 1612 1613 PetscFunctionBegin; 1614 its = its*lits; 1615 1616 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1617 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1618 a->fshift = fshift; 1619 a->omega = omega; 1620 1621 diag = a->diag; 1622 t = a->ssor_work; 1623 idiag = a->idiag; 1624 mdiag = a->mdiag; 1625 1626 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1627 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1628 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1629 if (flag == SOR_APPLY_UPPER) { 1630 /* apply (U + D/omega) to the vector */ 1631 bs = b; 1632 for (i=0; i<m; i++) { 1633 d = fshift + mdiag[i]; 1634 n = a->i[i+1] - diag[i] - 1; 1635 idx = a->j + diag[i] + 1; 1636 v = a->a + diag[i] + 1; 1637 sum = b[i]*d/omega; 1638 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1639 x[i] = sum; 1640 } 1641 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1642 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1643 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1644 PetscFunctionReturn(0); 1645 } 1646 1647 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1648 else if (flag & SOR_EISENSTAT) { 1649 /* Let A = L + U + D; where L is lower trianglar, 1650 U is upper triangular, E = D/omega; This routine applies 1651 1652 (L + E)^{-1} A (U + E)^{-1} 1653 1654 to a vector efficiently using Eisenstat's trick. 1655 */ 1656 scale = (2.0/omega) - 1.0; 1657 1658 /* x = (E + U)^{-1} b */ 1659 for (i=m-1; i>=0; i--) { 1660 n = a->i[i+1] - diag[i] - 1; 1661 idx = a->j + diag[i] + 1; 1662 v = a->a + diag[i] + 1; 1663 sum = b[i]; 1664 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1665 x[i] = sum*idiag[i]; 1666 } 1667 1668 /* t = b - (2*E - D)x */ 1669 v = a->a; 1670 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1671 1672 /* t = (E + L)^{-1}t */ 1673 ts = t; 1674 diag = a->diag; 1675 for (i=0; i<m; i++) { 1676 n = diag[i] - a->i[i]; 1677 idx = a->j + a->i[i]; 1678 v = a->a + a->i[i]; 1679 sum = t[i]; 1680 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1681 t[i] = sum*idiag[i]; 1682 /* x = x + t */ 1683 x[i] += t[i]; 1684 } 1685 1686 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1687 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1688 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1689 PetscFunctionReturn(0); 1690 } 1691 if (flag & SOR_ZERO_INITIAL_GUESS) { 1692 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1693 for (i=0; i<m; i++) { 1694 n = diag[i] - a->i[i]; 1695 idx = a->j + a->i[i]; 1696 v = a->a + a->i[i]; 1697 sum = b[i]; 1698 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1699 t[i] = sum; 1700 x[i] = sum*idiag[i]; 1701 } 1702 xb = t; 1703 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1704 } else xb = b; 1705 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1706 for (i=m-1; i>=0; i--) { 1707 n = a->i[i+1] - diag[i] - 1; 1708 idx = a->j + diag[i] + 1; 1709 v = a->a + diag[i] + 1; 1710 sum = xb[i]; 1711 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1712 if (xb == b) { 1713 x[i] = sum*idiag[i]; 1714 } else { 1715 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1716 } 1717 } 1718 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1719 } 1720 its--; 1721 } 1722 while (its--) { 1723 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1724 for (i=0; i<m; i++) { 1725 /* lower */ 1726 n = diag[i] - a->i[i]; 1727 idx = a->j + a->i[i]; 1728 v = a->a + a->i[i]; 1729 sum = b[i]; 1730 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1731 t[i] = sum; /* save application of the lower-triangular part */ 1732 /* upper */ 1733 n = a->i[i+1] - diag[i] - 1; 1734 idx = a->j + diag[i] + 1; 1735 v = a->a + diag[i] + 1; 1736 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1737 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1738 } 1739 xb = t; 1740 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1741 } else xb = b; 1742 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1743 for (i=m-1; i>=0; i--) { 1744 sum = xb[i]; 1745 if (xb == b) { 1746 /* whole matrix (no checkpointing available) */ 1747 n = a->i[i+1] - a->i[i]; 1748 idx = a->j + a->i[i]; 1749 v = a->a + a->i[i]; 1750 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1751 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1752 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1753 n = a->i[i+1] - diag[i] - 1; 1754 idx = a->j + diag[i] + 1; 1755 v = a->a + diag[i] + 1; 1756 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1757 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1758 } 1759 } 1760 if (xb == b) { 1761 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1762 } else { 1763 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1764 } 1765 } 1766 } 1767 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1768 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1769 PetscFunctionReturn(0); 1770 } 1771 1772 1773 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1774 { 1775 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1776 1777 PetscFunctionBegin; 1778 info->block_size = 1.0; 1779 info->nz_allocated = (double)a->maxnz; 1780 info->nz_used = (double)a->nz; 1781 info->nz_unneeded = (double)(a->maxnz - a->nz); 1782 info->assemblies = (double)A->num_ass; 1783 info->mallocs = (double)A->info.mallocs; 1784 info->memory = ((PetscObject)A)->mem; 1785 if (A->factortype) { 1786 info->fill_ratio_given = A->info.fill_ratio_given; 1787 info->fill_ratio_needed = A->info.fill_ratio_needed; 1788 info->factor_mallocs = A->info.factor_mallocs; 1789 } else { 1790 info->fill_ratio_given = 0; 1791 info->fill_ratio_needed = 0; 1792 info->factor_mallocs = 0; 1793 } 1794 PetscFunctionReturn(0); 1795 } 1796 1797 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1798 { 1799 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1800 PetscInt i,m = A->rmap->n - 1; 1801 PetscErrorCode ierr; 1802 const PetscScalar *xx; 1803 PetscScalar *bb; 1804 PetscInt d = 0; 1805 1806 PetscFunctionBegin; 1807 if (x && b) { 1808 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1809 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1810 for (i=0; i<N; i++) { 1811 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1812 bb[rows[i]] = diag*xx[rows[i]]; 1813 } 1814 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1815 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1816 } 1817 1818 if (a->keepnonzeropattern) { 1819 for (i=0; i<N; i++) { 1820 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1821 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1822 } 1823 if (diag != 0.0) { 1824 for (i=0; i<N; i++) { 1825 d = rows[i]; 1826 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); 1827 } 1828 for (i=0; i<N; i++) { 1829 a->a[a->diag[rows[i]]] = diag; 1830 } 1831 } 1832 } else { 1833 if (diag != 0.0) { 1834 for (i=0; i<N; i++) { 1835 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1836 if (a->ilen[rows[i]] > 0) { 1837 a->ilen[rows[i]] = 1; 1838 a->a[a->i[rows[i]]] = diag; 1839 a->j[a->i[rows[i]]] = rows[i]; 1840 } else { /* in case row was completely empty */ 1841 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1842 } 1843 } 1844 } else { 1845 for (i=0; i<N; i++) { 1846 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1847 a->ilen[rows[i]] = 0; 1848 } 1849 } 1850 A->nonzerostate++; 1851 } 1852 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1853 PetscFunctionReturn(0); 1854 } 1855 1856 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1857 { 1858 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1859 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1860 PetscErrorCode ierr; 1861 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1862 const PetscScalar *xx; 1863 PetscScalar *bb; 1864 1865 PetscFunctionBegin; 1866 if (x && b) { 1867 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1868 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1869 vecs = PETSC_TRUE; 1870 } 1871 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1872 for (i=0; i<N; i++) { 1873 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1874 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1875 1876 zeroed[rows[i]] = PETSC_TRUE; 1877 } 1878 for (i=0; i<A->rmap->n; i++) { 1879 if (!zeroed[i]) { 1880 for (j=a->i[i]; j<a->i[i+1]; j++) { 1881 if (zeroed[a->j[j]]) { 1882 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1883 a->a[j] = 0.0; 1884 } 1885 } 1886 } else if (vecs) bb[i] = diag*xx[i]; 1887 } 1888 if (x && b) { 1889 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1890 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1891 } 1892 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1893 if (diag != 0.0) { 1894 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1895 if (missing) { 1896 if (a->nonew) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1897 else { 1898 for (i=0; i<N; i++) { 1899 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1900 } 1901 } 1902 } else { 1903 for (i=0; i<N; i++) { 1904 a->a[a->diag[rows[i]]] = diag; 1905 } 1906 } 1907 } 1908 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1909 PetscFunctionReturn(0); 1910 } 1911 1912 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1913 { 1914 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1915 PetscInt *itmp; 1916 1917 PetscFunctionBegin; 1918 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1919 1920 *nz = a->i[row+1] - a->i[row]; 1921 if (v) *v = a->a + a->i[row]; 1922 if (idx) { 1923 itmp = a->j + a->i[row]; 1924 if (*nz) *idx = itmp; 1925 else *idx = 0; 1926 } 1927 PetscFunctionReturn(0); 1928 } 1929 1930 /* remove this function? */ 1931 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1932 { 1933 PetscFunctionBegin; 1934 PetscFunctionReturn(0); 1935 } 1936 1937 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1938 { 1939 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1940 MatScalar *v = a->a; 1941 PetscReal sum = 0.0; 1942 PetscErrorCode ierr; 1943 PetscInt i,j; 1944 1945 PetscFunctionBegin; 1946 if (type == NORM_FROBENIUS) { 1947 #if defined(PETSC_USE_REAL___FP16) 1948 PetscBLASInt one = 1,nz = a->nz; 1949 *nrm = BLASnrm2_(&nz,v,&one); 1950 #else 1951 for (i=0; i<a->nz; i++) { 1952 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1953 } 1954 *nrm = PetscSqrtReal(sum); 1955 #endif 1956 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 1957 } else if (type == NORM_1) { 1958 PetscReal *tmp; 1959 PetscInt *jj = a->j; 1960 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 1961 *nrm = 0.0; 1962 for (j=0; j<a->nz; j++) { 1963 tmp[*jj++] += PetscAbsScalar(*v); v++; 1964 } 1965 for (j=0; j<A->cmap->n; j++) { 1966 if (tmp[j] > *nrm) *nrm = tmp[j]; 1967 } 1968 ierr = PetscFree(tmp);CHKERRQ(ierr); 1969 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1970 } else if (type == NORM_INFINITY) { 1971 *nrm = 0.0; 1972 for (j=0; j<A->rmap->n; j++) { 1973 v = a->a + a->i[j]; 1974 sum = 0.0; 1975 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1976 sum += PetscAbsScalar(*v); v++; 1977 } 1978 if (sum > *nrm) *nrm = sum; 1979 } 1980 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1981 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1982 PetscFunctionReturn(0); 1983 } 1984 1985 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 1986 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 1987 { 1988 PetscErrorCode ierr; 1989 PetscInt i,j,anzj; 1990 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 1991 PetscInt an=A->cmap->N,am=A->rmap->N; 1992 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 1993 1994 PetscFunctionBegin; 1995 /* Allocate space for symbolic transpose info and work array */ 1996 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 1997 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 1998 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 1999 2000 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2001 /* Note: offset by 1 for fast conversion into csr format. */ 2002 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2003 /* Form ati for csr format of A^T. */ 2004 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2005 2006 /* Copy ati into atfill so we have locations of the next free space in atj */ 2007 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2008 2009 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2010 for (i=0;i<am;i++) { 2011 anzj = ai[i+1] - ai[i]; 2012 for (j=0;j<anzj;j++) { 2013 atj[atfill[*aj]] = i; 2014 atfill[*aj++] += 1; 2015 } 2016 } 2017 2018 /* Clean up temporary space and complete requests. */ 2019 ierr = PetscFree(atfill);CHKERRQ(ierr); 2020 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2021 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2022 2023 b = (Mat_SeqAIJ*)((*B)->data); 2024 b->free_a = PETSC_FALSE; 2025 b->free_ij = PETSC_TRUE; 2026 b->nonew = 0; 2027 PetscFunctionReturn(0); 2028 } 2029 2030 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2031 { 2032 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2033 Mat C; 2034 PetscErrorCode ierr; 2035 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2036 MatScalar *array = a->a; 2037 2038 PetscFunctionBegin; 2039 if (reuse == MAT_INPLACE_MATRIX && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 2040 2041 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 2042 ierr = PetscCalloc1(1+A->cmap->n,&col);CHKERRQ(ierr); 2043 2044 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2045 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2046 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2047 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2048 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2049 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2050 ierr = PetscFree(col);CHKERRQ(ierr); 2051 } else { 2052 C = *B; 2053 } 2054 2055 for (i=0; i<m; i++) { 2056 len = ai[i+1]-ai[i]; 2057 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2058 array += len; 2059 aj += len; 2060 } 2061 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2062 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2063 2064 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 2065 *B = C; 2066 } else { 2067 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 2068 } 2069 PetscFunctionReturn(0); 2070 } 2071 2072 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2073 { 2074 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2075 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2076 MatScalar *va,*vb; 2077 PetscErrorCode ierr; 2078 PetscInt ma,na,mb,nb, i; 2079 2080 PetscFunctionBegin; 2081 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2082 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2083 if (ma!=nb || na!=mb) { 2084 *f = PETSC_FALSE; 2085 PetscFunctionReturn(0); 2086 } 2087 aii = aij->i; bii = bij->i; 2088 adx = aij->j; bdx = bij->j; 2089 va = aij->a; vb = bij->a; 2090 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2091 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2092 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2093 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2094 2095 *f = PETSC_TRUE; 2096 for (i=0; i<ma; i++) { 2097 while (aptr[i]<aii[i+1]) { 2098 PetscInt idc,idr; 2099 PetscScalar vc,vr; 2100 /* column/row index/value */ 2101 idc = adx[aptr[i]]; 2102 idr = bdx[bptr[idc]]; 2103 vc = va[aptr[i]]; 2104 vr = vb[bptr[idc]]; 2105 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2106 *f = PETSC_FALSE; 2107 goto done; 2108 } else { 2109 aptr[i]++; 2110 if (B || i!=idc) bptr[idc]++; 2111 } 2112 } 2113 } 2114 done: 2115 ierr = PetscFree(aptr);CHKERRQ(ierr); 2116 ierr = PetscFree(bptr);CHKERRQ(ierr); 2117 PetscFunctionReturn(0); 2118 } 2119 2120 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2121 { 2122 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2123 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2124 MatScalar *va,*vb; 2125 PetscErrorCode ierr; 2126 PetscInt ma,na,mb,nb, i; 2127 2128 PetscFunctionBegin; 2129 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2130 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2131 if (ma!=nb || na!=mb) { 2132 *f = PETSC_FALSE; 2133 PetscFunctionReturn(0); 2134 } 2135 aii = aij->i; bii = bij->i; 2136 adx = aij->j; bdx = bij->j; 2137 va = aij->a; vb = bij->a; 2138 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2139 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2140 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2141 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2142 2143 *f = PETSC_TRUE; 2144 for (i=0; i<ma; i++) { 2145 while (aptr[i]<aii[i+1]) { 2146 PetscInt idc,idr; 2147 PetscScalar vc,vr; 2148 /* column/row index/value */ 2149 idc = adx[aptr[i]]; 2150 idr = bdx[bptr[idc]]; 2151 vc = va[aptr[i]]; 2152 vr = vb[bptr[idc]]; 2153 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2154 *f = PETSC_FALSE; 2155 goto done; 2156 } else { 2157 aptr[i]++; 2158 if (B || i!=idc) bptr[idc]++; 2159 } 2160 } 2161 } 2162 done: 2163 ierr = PetscFree(aptr);CHKERRQ(ierr); 2164 ierr = PetscFree(bptr);CHKERRQ(ierr); 2165 PetscFunctionReturn(0); 2166 } 2167 2168 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2169 { 2170 PetscErrorCode ierr; 2171 2172 PetscFunctionBegin; 2173 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2174 PetscFunctionReturn(0); 2175 } 2176 2177 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2178 { 2179 PetscErrorCode ierr; 2180 2181 PetscFunctionBegin; 2182 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2183 PetscFunctionReturn(0); 2184 } 2185 2186 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2187 { 2188 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2189 const PetscScalar *l,*r; 2190 PetscScalar x; 2191 MatScalar *v; 2192 PetscErrorCode ierr; 2193 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2194 const PetscInt *jj; 2195 2196 PetscFunctionBegin; 2197 if (ll) { 2198 /* The local size is used so that VecMPI can be passed to this routine 2199 by MatDiagonalScale_MPIAIJ */ 2200 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2201 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2202 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 2203 v = a->a; 2204 for (i=0; i<m; i++) { 2205 x = l[i]; 2206 M = a->i[i+1] - a->i[i]; 2207 for (j=0; j<M; j++) (*v++) *= x; 2208 } 2209 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 2210 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2211 } 2212 if (rr) { 2213 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2214 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2215 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 2216 v = a->a; jj = a->j; 2217 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2218 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 2219 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2220 } 2221 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2222 PetscFunctionReturn(0); 2223 } 2224 2225 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2226 { 2227 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2228 PetscErrorCode ierr; 2229 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2230 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2231 const PetscInt *irow,*icol; 2232 PetscInt nrows,ncols; 2233 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2234 MatScalar *a_new,*mat_a; 2235 Mat C; 2236 PetscBool stride; 2237 2238 PetscFunctionBegin; 2239 2240 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2241 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2242 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2243 2244 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2245 if (stride) { 2246 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2247 } else { 2248 first = 0; 2249 step = 0; 2250 } 2251 if (stride && step == 1) { 2252 /* special case of contiguous rows */ 2253 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2254 /* loop over new rows determining lens and starting points */ 2255 for (i=0; i<nrows; i++) { 2256 kstart = ai[irow[i]]; 2257 kend = kstart + ailen[irow[i]]; 2258 starts[i] = kstart; 2259 for (k=kstart; k<kend; k++) { 2260 if (aj[k] >= first) { 2261 starts[i] = k; 2262 break; 2263 } 2264 } 2265 sum = 0; 2266 while (k < kend) { 2267 if (aj[k++] >= first+ncols) break; 2268 sum++; 2269 } 2270 lens[i] = sum; 2271 } 2272 /* create submatrix */ 2273 if (scall == MAT_REUSE_MATRIX) { 2274 PetscInt n_cols,n_rows; 2275 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2276 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2277 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2278 C = *B; 2279 } else { 2280 PetscInt rbs,cbs; 2281 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2282 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2283 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2284 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2285 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2286 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2287 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2288 } 2289 c = (Mat_SeqAIJ*)C->data; 2290 2291 /* loop over rows inserting into submatrix */ 2292 a_new = c->a; 2293 j_new = c->j; 2294 i_new = c->i; 2295 2296 for (i=0; i<nrows; i++) { 2297 ii = starts[i]; 2298 lensi = lens[i]; 2299 for (k=0; k<lensi; k++) { 2300 *j_new++ = aj[ii+k] - first; 2301 } 2302 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2303 a_new += lensi; 2304 i_new[i+1] = i_new[i] + lensi; 2305 c->ilen[i] = lensi; 2306 } 2307 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2308 } else { 2309 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2310 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2311 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2312 for (i=0; i<ncols; i++) { 2313 #if defined(PETSC_USE_DEBUG) 2314 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); 2315 #endif 2316 smap[icol[i]] = i+1; 2317 } 2318 2319 /* determine lens of each row */ 2320 for (i=0; i<nrows; i++) { 2321 kstart = ai[irow[i]]; 2322 kend = kstart + a->ilen[irow[i]]; 2323 lens[i] = 0; 2324 for (k=kstart; k<kend; k++) { 2325 if (smap[aj[k]]) { 2326 lens[i]++; 2327 } 2328 } 2329 } 2330 /* Create and fill new matrix */ 2331 if (scall == MAT_REUSE_MATRIX) { 2332 PetscBool equal; 2333 2334 c = (Mat_SeqAIJ*)((*B)->data); 2335 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2336 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2337 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2338 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2339 C = *B; 2340 } else { 2341 PetscInt rbs,cbs; 2342 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2343 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2344 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2345 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2346 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2347 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2348 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2349 } 2350 c = (Mat_SeqAIJ*)(C->data); 2351 for (i=0; i<nrows; i++) { 2352 row = irow[i]; 2353 kstart = ai[row]; 2354 kend = kstart + a->ilen[row]; 2355 mat_i = c->i[i]; 2356 mat_j = c->j + mat_i; 2357 mat_a = c->a + mat_i; 2358 mat_ilen = c->ilen + i; 2359 for (k=kstart; k<kend; k++) { 2360 if ((tcol=smap[a->j[k]])) { 2361 *mat_j++ = tcol - 1; 2362 *mat_a++ = a->a[k]; 2363 (*mat_ilen)++; 2364 2365 } 2366 } 2367 } 2368 /* Free work space */ 2369 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2370 ierr = PetscFree(smap);CHKERRQ(ierr); 2371 ierr = PetscFree(lens);CHKERRQ(ierr); 2372 /* sort */ 2373 for (i = 0; i < nrows; i++) { 2374 PetscInt ilen; 2375 2376 mat_i = c->i[i]; 2377 mat_j = c->j + mat_i; 2378 mat_a = c->a + mat_i; 2379 ilen = c->ilen[i]; 2380 ierr = PetscSortIntWithScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2381 } 2382 } 2383 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2384 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2385 2386 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2387 *B = C; 2388 PetscFunctionReturn(0); 2389 } 2390 2391 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2392 { 2393 PetscErrorCode ierr; 2394 Mat B; 2395 2396 PetscFunctionBegin; 2397 if (scall == MAT_INITIAL_MATRIX) { 2398 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2399 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2400 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2401 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2402 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2403 *subMat = B; 2404 } else { 2405 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2406 } 2407 PetscFunctionReturn(0); 2408 } 2409 2410 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2411 { 2412 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2413 PetscErrorCode ierr; 2414 Mat outA; 2415 PetscBool row_identity,col_identity; 2416 2417 PetscFunctionBegin; 2418 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2419 2420 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2421 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2422 2423 outA = inA; 2424 outA->factortype = MAT_FACTOR_LU; 2425 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2426 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2427 2428 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2429 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2430 2431 a->row = row; 2432 2433 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2434 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2435 2436 a->col = col; 2437 2438 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2439 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2440 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2441 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2442 2443 if (!a->solve_work) { /* this matrix may have been factored before */ 2444 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2445 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2446 } 2447 2448 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2449 if (row_identity && col_identity) { 2450 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2451 } else { 2452 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2453 } 2454 PetscFunctionReturn(0); 2455 } 2456 2457 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2458 { 2459 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2460 PetscScalar oalpha = alpha; 2461 PetscErrorCode ierr; 2462 PetscBLASInt one = 1,bnz; 2463 2464 PetscFunctionBegin; 2465 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2466 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2467 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2468 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2469 PetscFunctionReturn(0); 2470 } 2471 2472 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2473 { 2474 PetscErrorCode ierr; 2475 PetscInt i; 2476 2477 PetscFunctionBegin; 2478 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2479 ierr = PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);CHKERRQ(ierr); 2480 2481 for (i=0; i<submatj->nrqr; ++i) { 2482 ierr = PetscFree(submatj->sbuf2[i]);CHKERRQ(ierr); 2483 } 2484 ierr = PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);CHKERRQ(ierr); 2485 2486 if (submatj->rbuf1) { 2487 ierr = PetscFree(submatj->rbuf1[0]);CHKERRQ(ierr); 2488 ierr = PetscFree(submatj->rbuf1);CHKERRQ(ierr); 2489 } 2490 2491 for (i=0; i<submatj->nrqs; ++i) { 2492 ierr = PetscFree(submatj->rbuf3[i]);CHKERRQ(ierr); 2493 } 2494 ierr = PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);CHKERRQ(ierr); 2495 ierr = PetscFree(submatj->pa);CHKERRQ(ierr); 2496 } 2497 2498 #if defined(PETSC_USE_CTABLE) 2499 ierr = PetscTableDestroy((PetscTable*)&submatj->rmap);CHKERRQ(ierr); 2500 if (submatj->cmap_loc) {ierr = PetscFree(submatj->cmap_loc);CHKERRQ(ierr);} 2501 ierr = PetscFree(submatj->rmap_loc);CHKERRQ(ierr); 2502 #else 2503 ierr = PetscFree(submatj->rmap);CHKERRQ(ierr); 2504 #endif 2505 2506 if (!submatj->allcolumns) { 2507 #if defined(PETSC_USE_CTABLE) 2508 ierr = PetscTableDestroy((PetscTable*)&submatj->cmap);CHKERRQ(ierr); 2509 #else 2510 ierr = PetscFree(submatj->cmap);CHKERRQ(ierr); 2511 #endif 2512 } 2513 ierr = PetscFree(submatj->row2proc);CHKERRQ(ierr); 2514 2515 ierr = PetscFree(submatj);CHKERRQ(ierr); 2516 PetscFunctionReturn(0); 2517 } 2518 2519 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2520 { 2521 PetscErrorCode ierr; 2522 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2523 Mat_SubSppt *submatj = c->submatis1; 2524 2525 PetscFunctionBegin; 2526 ierr = submatj->destroy(C);CHKERRQ(ierr); 2527 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2528 PetscFunctionReturn(0); 2529 } 2530 2531 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2532 { 2533 PetscErrorCode ierr; 2534 PetscInt i; 2535 Mat C; 2536 Mat_SeqAIJ *c; 2537 Mat_SubSppt *submatj; 2538 2539 PetscFunctionBegin; 2540 for (i=0; i<n; i++) { 2541 C = (*mat)[i]; 2542 c = (Mat_SeqAIJ*)C->data; 2543 submatj = c->submatis1; 2544 if (submatj) { 2545 if (--((PetscObject)C)->refct <= 0) { 2546 ierr = (submatj->destroy)(C);CHKERRQ(ierr); 2547 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2548 ierr = PetscLayoutDestroy(&C->rmap);CHKERRQ(ierr); 2549 ierr = PetscLayoutDestroy(&C->cmap);CHKERRQ(ierr); 2550 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 2551 } 2552 } else { 2553 ierr = MatDestroy(&C);CHKERRQ(ierr); 2554 } 2555 } 2556 2557 ierr = PetscFree(*mat);CHKERRQ(ierr); 2558 PetscFunctionReturn(0); 2559 } 2560 2561 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2562 { 2563 PetscErrorCode ierr; 2564 PetscInt i; 2565 2566 PetscFunctionBegin; 2567 if (scall == MAT_INITIAL_MATRIX) { 2568 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2569 } 2570 2571 for (i=0; i<n; i++) { 2572 ierr = MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2573 } 2574 PetscFunctionReturn(0); 2575 } 2576 2577 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2578 { 2579 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2580 PetscErrorCode ierr; 2581 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2582 const PetscInt *idx; 2583 PetscInt start,end,*ai,*aj; 2584 PetscBT table; 2585 2586 PetscFunctionBegin; 2587 m = A->rmap->n; 2588 ai = a->i; 2589 aj = a->j; 2590 2591 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2592 2593 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2594 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2595 2596 for (i=0; i<is_max; i++) { 2597 /* Initialize the two local arrays */ 2598 isz = 0; 2599 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2600 2601 /* Extract the indices, assume there can be duplicate entries */ 2602 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2603 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2604 2605 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2606 for (j=0; j<n; ++j) { 2607 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2608 } 2609 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2610 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2611 2612 k = 0; 2613 for (j=0; j<ov; j++) { /* for each overlap */ 2614 n = isz; 2615 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2616 row = nidx[k]; 2617 start = ai[row]; 2618 end = ai[row+1]; 2619 for (l = start; l<end; l++) { 2620 val = aj[l]; 2621 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2622 } 2623 } 2624 } 2625 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2626 } 2627 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2628 ierr = PetscFree(nidx);CHKERRQ(ierr); 2629 PetscFunctionReturn(0); 2630 } 2631 2632 /* -------------------------------------------------------------- */ 2633 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2634 { 2635 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2636 PetscErrorCode ierr; 2637 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2638 const PetscInt *row,*col; 2639 PetscInt *cnew,j,*lens; 2640 IS icolp,irowp; 2641 PetscInt *cwork = NULL; 2642 PetscScalar *vwork = NULL; 2643 2644 PetscFunctionBegin; 2645 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2646 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2647 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2648 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2649 2650 /* determine lengths of permuted rows */ 2651 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2652 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2653 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2654 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2655 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2656 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2657 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2658 ierr = PetscFree(lens);CHKERRQ(ierr); 2659 2660 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2661 for (i=0; i<m; i++) { 2662 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2663 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2664 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2665 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2666 } 2667 ierr = PetscFree(cnew);CHKERRQ(ierr); 2668 2669 (*B)->assembled = PETSC_FALSE; 2670 2671 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2672 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2673 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2674 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2675 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2676 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2677 PetscFunctionReturn(0); 2678 } 2679 2680 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2681 { 2682 PetscErrorCode ierr; 2683 2684 PetscFunctionBegin; 2685 /* If the two matrices have the same copy implementation, use fast copy. */ 2686 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2687 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2688 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2689 2690 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"); 2691 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2692 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2693 } else { 2694 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2695 } 2696 PetscFunctionReturn(0); 2697 } 2698 2699 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2700 { 2701 PetscErrorCode ierr; 2702 2703 PetscFunctionBegin; 2704 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2705 PetscFunctionReturn(0); 2706 } 2707 2708 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2709 { 2710 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2711 2712 PetscFunctionBegin; 2713 *array = a->a; 2714 PetscFunctionReturn(0); 2715 } 2716 2717 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2718 { 2719 PetscFunctionBegin; 2720 PetscFunctionReturn(0); 2721 } 2722 2723 /* 2724 Computes the number of nonzeros per row needed for preallocation when X and Y 2725 have different nonzero structure. 2726 */ 2727 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2728 { 2729 PetscInt i,j,k,nzx,nzy; 2730 2731 PetscFunctionBegin; 2732 /* Set the number of nonzeros in the new matrix */ 2733 for (i=0; i<m; i++) { 2734 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2735 nzx = xi[i+1] - xi[i]; 2736 nzy = yi[i+1] - yi[i]; 2737 nnz[i] = 0; 2738 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2739 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2740 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2741 nnz[i]++; 2742 } 2743 for (; k<nzy; k++) nnz[i]++; 2744 } 2745 PetscFunctionReturn(0); 2746 } 2747 2748 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2749 { 2750 PetscInt m = Y->rmap->N; 2751 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2752 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2753 PetscErrorCode ierr; 2754 2755 PetscFunctionBegin; 2756 /* Set the number of nonzeros in the new matrix */ 2757 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2758 PetscFunctionReturn(0); 2759 } 2760 2761 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2762 { 2763 PetscErrorCode ierr; 2764 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2765 PetscBLASInt one=1,bnz; 2766 2767 PetscFunctionBegin; 2768 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2769 if (str == SAME_NONZERO_PATTERN) { 2770 PetscScalar alpha = a; 2771 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2772 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2773 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2774 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2775 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2776 } else { 2777 Mat B; 2778 PetscInt *nnz; 2779 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2780 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2781 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2782 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2783 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2784 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2785 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2786 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2787 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2788 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2789 ierr = PetscFree(nnz);CHKERRQ(ierr); 2790 } 2791 PetscFunctionReturn(0); 2792 } 2793 2794 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2795 { 2796 #if defined(PETSC_USE_COMPLEX) 2797 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2798 PetscInt i,nz; 2799 PetscScalar *a; 2800 2801 PetscFunctionBegin; 2802 nz = aij->nz; 2803 a = aij->a; 2804 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2805 #else 2806 PetscFunctionBegin; 2807 #endif 2808 PetscFunctionReturn(0); 2809 } 2810 2811 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2812 { 2813 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2814 PetscErrorCode ierr; 2815 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2816 PetscReal atmp; 2817 PetscScalar *x; 2818 MatScalar *aa; 2819 2820 PetscFunctionBegin; 2821 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2822 aa = a->a; 2823 ai = a->i; 2824 aj = a->j; 2825 2826 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2827 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2828 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2829 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2830 for (i=0; i<m; i++) { 2831 ncols = ai[1] - ai[0]; ai++; 2832 x[i] = 0.0; 2833 for (j=0; j<ncols; j++) { 2834 atmp = PetscAbsScalar(*aa); 2835 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2836 aa++; aj++; 2837 } 2838 } 2839 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2840 PetscFunctionReturn(0); 2841 } 2842 2843 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2844 { 2845 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2846 PetscErrorCode ierr; 2847 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2848 PetscScalar *x; 2849 MatScalar *aa; 2850 2851 PetscFunctionBegin; 2852 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2853 aa = a->a; 2854 ai = a->i; 2855 aj = a->j; 2856 2857 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2858 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2859 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2860 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2861 for (i=0; i<m; i++) { 2862 ncols = ai[1] - ai[0]; ai++; 2863 if (ncols == A->cmap->n) { /* row is dense */ 2864 x[i] = *aa; if (idx) idx[i] = 0; 2865 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2866 x[i] = 0.0; 2867 if (idx) { 2868 idx[i] = 0; /* in case ncols is zero */ 2869 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2870 if (aj[j] > j) { 2871 idx[i] = j; 2872 break; 2873 } 2874 } 2875 } 2876 } 2877 for (j=0; j<ncols; j++) { 2878 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2879 aa++; aj++; 2880 } 2881 } 2882 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2883 PetscFunctionReturn(0); 2884 } 2885 2886 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2887 { 2888 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2889 PetscErrorCode ierr; 2890 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2891 PetscReal atmp; 2892 PetscScalar *x; 2893 MatScalar *aa; 2894 2895 PetscFunctionBegin; 2896 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2897 aa = a->a; 2898 ai = a->i; 2899 aj = a->j; 2900 2901 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2902 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2903 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2904 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); 2905 for (i=0; i<m; i++) { 2906 ncols = ai[1] - ai[0]; ai++; 2907 if (ncols) { 2908 /* Get first nonzero */ 2909 for (j = 0; j < ncols; j++) { 2910 atmp = PetscAbsScalar(aa[j]); 2911 if (atmp > 1.0e-12) { 2912 x[i] = atmp; 2913 if (idx) idx[i] = aj[j]; 2914 break; 2915 } 2916 } 2917 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2918 } else { 2919 x[i] = 0.0; if (idx) idx[i] = 0; 2920 } 2921 for (j = 0; j < ncols; j++) { 2922 atmp = PetscAbsScalar(*aa); 2923 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2924 aa++; aj++; 2925 } 2926 } 2927 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2928 PetscFunctionReturn(0); 2929 } 2930 2931 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2932 { 2933 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2934 PetscErrorCode ierr; 2935 PetscInt i,j,m = A->rmap->n,ncols,n; 2936 const PetscInt *ai,*aj; 2937 PetscScalar *x; 2938 const MatScalar *aa; 2939 2940 PetscFunctionBegin; 2941 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2942 aa = a->a; 2943 ai = a->i; 2944 aj = a->j; 2945 2946 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2947 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2948 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2949 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2950 for (i=0; i<m; i++) { 2951 ncols = ai[1] - ai[0]; ai++; 2952 if (ncols == A->cmap->n) { /* row is dense */ 2953 x[i] = *aa; if (idx) idx[i] = 0; 2954 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2955 x[i] = 0.0; 2956 if (idx) { /* find first implicit 0.0 in the row */ 2957 idx[i] = 0; /* in case ncols is zero */ 2958 for (j=0; j<ncols; j++) { 2959 if (aj[j] > j) { 2960 idx[i] = j; 2961 break; 2962 } 2963 } 2964 } 2965 } 2966 for (j=0; j<ncols; j++) { 2967 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2968 aa++; aj++; 2969 } 2970 } 2971 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2972 PetscFunctionReturn(0); 2973 } 2974 2975 #include <petscblaslapack.h> 2976 #include <petsc/private/kernels/blockinvert.h> 2977 2978 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 2979 { 2980 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2981 PetscErrorCode ierr; 2982 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2983 MatScalar *diag,work[25],*v_work; 2984 PetscReal shift = 0.0; 2985 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 2986 2987 PetscFunctionBegin; 2988 allowzeropivot = PetscNot(A->erroriffailure); 2989 if (a->ibdiagvalid) { 2990 if (values) *values = a->ibdiag; 2991 PetscFunctionReturn(0); 2992 } 2993 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2994 if (!a->ibdiag) { 2995 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 2996 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2997 } 2998 diag = a->ibdiag; 2999 if (values) *values = a->ibdiag; 3000 /* factor and invert each block */ 3001 switch (bs) { 3002 case 1: 3003 for (i=0; i<mbs; i++) { 3004 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3005 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3006 if (allowzeropivot) { 3007 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3008 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3009 A->factorerror_zeropivot_row = i; 3010 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 3011 } 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); 3012 } 3013 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3014 } 3015 break; 3016 case 2: 3017 for (i=0; i<mbs; i++) { 3018 ij[0] = 2*i; ij[1] = 2*i + 1; 3019 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3020 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3021 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3022 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3023 diag += 4; 3024 } 3025 break; 3026 case 3: 3027 for (i=0; i<mbs; i++) { 3028 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3029 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3030 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3031 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3032 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3033 diag += 9; 3034 } 3035 break; 3036 case 4: 3037 for (i=0; i<mbs; i++) { 3038 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3039 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3040 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3041 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3042 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3043 diag += 16; 3044 } 3045 break; 3046 case 5: 3047 for (i=0; i<mbs; i++) { 3048 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3049 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3050 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3051 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3052 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3053 diag += 25; 3054 } 3055 break; 3056 case 6: 3057 for (i=0; i<mbs; i++) { 3058 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; 3059 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3060 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3061 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3062 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3063 diag += 36; 3064 } 3065 break; 3066 case 7: 3067 for (i=0; i<mbs; i++) { 3068 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; 3069 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3070 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3071 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3072 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3073 diag += 49; 3074 } 3075 break; 3076 default: 3077 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3078 for (i=0; i<mbs; i++) { 3079 for (j=0; j<bs; j++) { 3080 IJ[j] = bs*i + j; 3081 } 3082 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3083 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3084 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3085 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3086 diag += bs2; 3087 } 3088 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3089 } 3090 a->ibdiagvalid = PETSC_TRUE; 3091 PetscFunctionReturn(0); 3092 } 3093 3094 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3095 { 3096 PetscErrorCode ierr; 3097 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3098 PetscScalar a; 3099 PetscInt m,n,i,j,col; 3100 3101 PetscFunctionBegin; 3102 if (!x->assembled) { 3103 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3104 for (i=0; i<m; i++) { 3105 for (j=0; j<aij->imax[i]; j++) { 3106 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3107 col = (PetscInt)(n*PetscRealPart(a)); 3108 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3109 } 3110 } 3111 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3112 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3113 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3114 PetscFunctionReturn(0); 3115 } 3116 3117 PetscErrorCode MatShift_SeqAIJ(Mat Y,PetscScalar a) 3118 { 3119 PetscErrorCode ierr; 3120 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)Y->data; 3121 3122 PetscFunctionBegin; 3123 if (!Y->preallocated || !aij->nz) { 3124 ierr = MatSeqAIJSetPreallocation(Y,1,NULL);CHKERRQ(ierr); 3125 } 3126 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 3127 PetscFunctionReturn(0); 3128 } 3129 3130 /* -------------------------------------------------------------------*/ 3131 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3132 MatGetRow_SeqAIJ, 3133 MatRestoreRow_SeqAIJ, 3134 MatMult_SeqAIJ, 3135 /* 4*/ MatMultAdd_SeqAIJ, 3136 MatMultTranspose_SeqAIJ, 3137 MatMultTransposeAdd_SeqAIJ, 3138 0, 3139 0, 3140 0, 3141 /* 10*/ 0, 3142 MatLUFactor_SeqAIJ, 3143 0, 3144 MatSOR_SeqAIJ, 3145 MatTranspose_SeqAIJ, 3146 /*1 5*/ MatGetInfo_SeqAIJ, 3147 MatEqual_SeqAIJ, 3148 MatGetDiagonal_SeqAIJ, 3149 MatDiagonalScale_SeqAIJ, 3150 MatNorm_SeqAIJ, 3151 /* 20*/ 0, 3152 MatAssemblyEnd_SeqAIJ, 3153 MatSetOption_SeqAIJ, 3154 MatZeroEntries_SeqAIJ, 3155 /* 24*/ MatZeroRows_SeqAIJ, 3156 0, 3157 0, 3158 0, 3159 0, 3160 /* 29*/ MatSetUp_SeqAIJ, 3161 0, 3162 0, 3163 0, 3164 0, 3165 /* 34*/ MatDuplicate_SeqAIJ, 3166 0, 3167 0, 3168 MatILUFactor_SeqAIJ, 3169 0, 3170 /* 39*/ MatAXPY_SeqAIJ, 3171 MatCreateSubMatrices_SeqAIJ, 3172 MatIncreaseOverlap_SeqAIJ, 3173 MatGetValues_SeqAIJ, 3174 MatCopy_SeqAIJ, 3175 /* 44*/ MatGetRowMax_SeqAIJ, 3176 MatScale_SeqAIJ, 3177 MatShift_SeqAIJ, 3178 MatDiagonalSet_SeqAIJ, 3179 MatZeroRowsColumns_SeqAIJ, 3180 /* 49*/ MatSetRandom_SeqAIJ, 3181 MatGetRowIJ_SeqAIJ, 3182 MatRestoreRowIJ_SeqAIJ, 3183 MatGetColumnIJ_SeqAIJ, 3184 MatRestoreColumnIJ_SeqAIJ, 3185 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3186 0, 3187 0, 3188 MatPermute_SeqAIJ, 3189 0, 3190 /* 59*/ 0, 3191 MatDestroy_SeqAIJ, 3192 MatView_SeqAIJ, 3193 0, 3194 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3195 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3196 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3197 0, 3198 0, 3199 0, 3200 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3201 MatGetRowMinAbs_SeqAIJ, 3202 0, 3203 0, 3204 0, 3205 /* 74*/ 0, 3206 MatFDColoringApply_AIJ, 3207 0, 3208 0, 3209 0, 3210 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3211 0, 3212 0, 3213 0, 3214 MatLoad_SeqAIJ, 3215 /* 84*/ MatIsSymmetric_SeqAIJ, 3216 MatIsHermitian_SeqAIJ, 3217 0, 3218 0, 3219 0, 3220 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3221 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3222 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3223 MatPtAP_SeqAIJ_SeqAIJ, 3224 MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy, 3225 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3226 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3227 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3228 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3229 0, 3230 /* 99*/ 0, 3231 0, 3232 0, 3233 MatConjugate_SeqAIJ, 3234 0, 3235 /*104*/ MatSetValuesRow_SeqAIJ, 3236 MatRealPart_SeqAIJ, 3237 MatImaginaryPart_SeqAIJ, 3238 0, 3239 0, 3240 /*109*/ MatMatSolve_SeqAIJ, 3241 0, 3242 MatGetRowMin_SeqAIJ, 3243 0, 3244 MatMissingDiagonal_SeqAIJ, 3245 /*114*/ 0, 3246 0, 3247 0, 3248 0, 3249 0, 3250 /*119*/ 0, 3251 0, 3252 0, 3253 0, 3254 MatGetMultiProcBlock_SeqAIJ, 3255 /*124*/ MatFindNonzeroRows_SeqAIJ, 3256 MatGetColumnNorms_SeqAIJ, 3257 MatInvertBlockDiagonal_SeqAIJ, 3258 0, 3259 0, 3260 /*129*/ 0, 3261 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3262 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3263 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3264 MatTransposeColoringCreate_SeqAIJ, 3265 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3266 MatTransColoringApplyDenToSp_SeqAIJ, 3267 MatRARt_SeqAIJ_SeqAIJ, 3268 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3269 MatRARtNumeric_SeqAIJ_SeqAIJ, 3270 /*139*/0, 3271 0, 3272 0, 3273 MatFDColoringSetUp_SeqXAIJ, 3274 MatFindOffBlockDiagonalEntries_SeqAIJ, 3275 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3276 MatDestroySubMatrices_SeqAIJ 3277 }; 3278 3279 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3280 { 3281 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3282 PetscInt i,nz,n; 3283 3284 PetscFunctionBegin; 3285 nz = aij->maxnz; 3286 n = mat->rmap->n; 3287 for (i=0; i<nz; i++) { 3288 aij->j[i] = indices[i]; 3289 } 3290 aij->nz = nz; 3291 for (i=0; i<n; i++) { 3292 aij->ilen[i] = aij->imax[i]; 3293 } 3294 PetscFunctionReturn(0); 3295 } 3296 3297 /*@ 3298 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3299 in the matrix. 3300 3301 Input Parameters: 3302 + mat - the SeqAIJ matrix 3303 - indices - the column indices 3304 3305 Level: advanced 3306 3307 Notes: 3308 This can be called if you have precomputed the nonzero structure of the 3309 matrix and want to provide it to the matrix object to improve the performance 3310 of the MatSetValues() operation. 3311 3312 You MUST have set the correct numbers of nonzeros per row in the call to 3313 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3314 3315 MUST be called before any calls to MatSetValues(); 3316 3317 The indices should start with zero, not one. 3318 3319 @*/ 3320 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3321 { 3322 PetscErrorCode ierr; 3323 3324 PetscFunctionBegin; 3325 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3326 PetscValidPointer(indices,2); 3327 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3328 PetscFunctionReturn(0); 3329 } 3330 3331 /* ----------------------------------------------------------------------------------------*/ 3332 3333 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3334 { 3335 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3336 PetscErrorCode ierr; 3337 size_t nz = aij->i[mat->rmap->n]; 3338 3339 PetscFunctionBegin; 3340 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3341 3342 /* allocate space for values if not already there */ 3343 if (!aij->saved_values) { 3344 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3345 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3346 } 3347 3348 /* copy values over */ 3349 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3350 PetscFunctionReturn(0); 3351 } 3352 3353 /*@ 3354 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3355 example, reuse of the linear part of a Jacobian, while recomputing the 3356 nonlinear portion. 3357 3358 Collect on Mat 3359 3360 Input Parameters: 3361 . mat - the matrix (currently only AIJ matrices support this option) 3362 3363 Level: advanced 3364 3365 Common Usage, with SNESSolve(): 3366 $ Create Jacobian matrix 3367 $ Set linear terms into matrix 3368 $ Apply boundary conditions to matrix, at this time matrix must have 3369 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3370 $ boundary conditions again will not change the nonzero structure 3371 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3372 $ ierr = MatStoreValues(mat); 3373 $ Call SNESSetJacobian() with matrix 3374 $ In your Jacobian routine 3375 $ ierr = MatRetrieveValues(mat); 3376 $ Set nonlinear terms in matrix 3377 3378 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3379 $ // build linear portion of Jacobian 3380 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3381 $ ierr = MatStoreValues(mat); 3382 $ loop over nonlinear iterations 3383 $ ierr = MatRetrieveValues(mat); 3384 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3385 $ // call MatAssemblyBegin/End() on matrix 3386 $ Solve linear system with Jacobian 3387 $ endloop 3388 3389 Notes: 3390 Matrix must already be assemblied before calling this routine 3391 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3392 calling this routine. 3393 3394 When this is called multiple times it overwrites the previous set of stored values 3395 and does not allocated additional space. 3396 3397 .seealso: MatRetrieveValues() 3398 3399 @*/ 3400 PetscErrorCode MatStoreValues(Mat mat) 3401 { 3402 PetscErrorCode ierr; 3403 3404 PetscFunctionBegin; 3405 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3406 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3407 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3408 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3409 PetscFunctionReturn(0); 3410 } 3411 3412 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3413 { 3414 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3415 PetscErrorCode ierr; 3416 PetscInt nz = aij->i[mat->rmap->n]; 3417 3418 PetscFunctionBegin; 3419 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3420 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3421 /* copy values over */ 3422 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3423 PetscFunctionReturn(0); 3424 } 3425 3426 /*@ 3427 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3428 example, reuse of the linear part of a Jacobian, while recomputing the 3429 nonlinear portion. 3430 3431 Collect on Mat 3432 3433 Input Parameters: 3434 . mat - the matrix (currently only AIJ matrices support this option) 3435 3436 Level: advanced 3437 3438 .seealso: MatStoreValues() 3439 3440 @*/ 3441 PetscErrorCode MatRetrieveValues(Mat mat) 3442 { 3443 PetscErrorCode ierr; 3444 3445 PetscFunctionBegin; 3446 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3447 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3448 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3449 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3450 PetscFunctionReturn(0); 3451 } 3452 3453 3454 /* --------------------------------------------------------------------------------*/ 3455 /*@C 3456 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3457 (the default parallel PETSc format). For good matrix assembly performance 3458 the user should preallocate the matrix storage by setting the parameter nz 3459 (or the array nnz). By setting these parameters accurately, performance 3460 during matrix assembly can be increased by more than a factor of 50. 3461 3462 Collective on MPI_Comm 3463 3464 Input Parameters: 3465 + comm - MPI communicator, set to PETSC_COMM_SELF 3466 . m - number of rows 3467 . n - number of columns 3468 . nz - number of nonzeros per row (same for all rows) 3469 - nnz - array containing the number of nonzeros in the various rows 3470 (possibly different for each row) or NULL 3471 3472 Output Parameter: 3473 . A - the matrix 3474 3475 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3476 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3477 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3478 3479 Notes: 3480 If nnz is given then nz is ignored 3481 3482 The AIJ format (also called the Yale sparse matrix format or 3483 compressed row storage), is fully compatible with standard Fortran 77 3484 storage. That is, the stored row and column indices can begin at 3485 either one (as in Fortran) or zero. See the users' manual for details. 3486 3487 Specify the preallocated storage with either nz or nnz (not both). 3488 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3489 allocation. For large problems you MUST preallocate memory or you 3490 will get TERRIBLE performance, see the users' manual chapter on matrices. 3491 3492 By default, this format uses inodes (identical nodes) when possible, to 3493 improve numerical efficiency of matrix-vector products and solves. We 3494 search for consecutive rows with the same nonzero structure, thereby 3495 reusing matrix information to achieve increased efficiency. 3496 3497 Options Database Keys: 3498 + -mat_no_inode - Do not use inodes 3499 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3500 3501 Level: intermediate 3502 3503 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3504 3505 @*/ 3506 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3507 { 3508 PetscErrorCode ierr; 3509 3510 PetscFunctionBegin; 3511 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3512 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3513 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3514 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3515 PetscFunctionReturn(0); 3516 } 3517 3518 /*@C 3519 MatSeqAIJSetPreallocation - For good matrix assembly performance 3520 the user should preallocate the matrix storage by setting the parameter nz 3521 (or the array nnz). By setting these parameters accurately, performance 3522 during matrix assembly can be increased by more than a factor of 50. 3523 3524 Collective on MPI_Comm 3525 3526 Input Parameters: 3527 + B - The matrix 3528 . nz - number of nonzeros per row (same for all rows) 3529 - nnz - array containing the number of nonzeros in the various rows 3530 (possibly different for each row) or NULL 3531 3532 Notes: 3533 If nnz is given then nz is ignored 3534 3535 The AIJ format (also called the Yale sparse matrix format or 3536 compressed row storage), is fully compatible with standard Fortran 77 3537 storage. That is, the stored row and column indices can begin at 3538 either one (as in Fortran) or zero. See the users' manual for details. 3539 3540 Specify the preallocated storage with either nz or nnz (not both). 3541 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3542 allocation. For large problems you MUST preallocate memory or you 3543 will get TERRIBLE performance, see the users' manual chapter on matrices. 3544 3545 You can call MatGetInfo() to get information on how effective the preallocation was; 3546 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3547 You can also run with the option -info and look for messages with the string 3548 malloc in them to see if additional memory allocation was needed. 3549 3550 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3551 entries or columns indices 3552 3553 By default, this format uses inodes (identical nodes) when possible, to 3554 improve numerical efficiency of matrix-vector products and solves. We 3555 search for consecutive rows with the same nonzero structure, thereby 3556 reusing matrix information to achieve increased efficiency. 3557 3558 Options Database Keys: 3559 + -mat_no_inode - Do not use inodes 3560 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3561 - -mat_aij_oneindex - Internally use indexing starting at 1 3562 rather than 0. Note that when calling MatSetValues(), 3563 the user still MUST index entries starting at 0! 3564 3565 Level: intermediate 3566 3567 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3568 3569 @*/ 3570 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3571 { 3572 PetscErrorCode ierr; 3573 3574 PetscFunctionBegin; 3575 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3576 PetscValidType(B,1); 3577 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3578 PetscFunctionReturn(0); 3579 } 3580 3581 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3582 { 3583 Mat_SeqAIJ *b; 3584 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3585 PetscErrorCode ierr; 3586 PetscInt i; 3587 3588 PetscFunctionBegin; 3589 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3590 if (nz == MAT_SKIP_ALLOCATION) { 3591 skipallocation = PETSC_TRUE; 3592 nz = 0; 3593 } 3594 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3595 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3596 3597 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3598 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3599 if (nnz) { 3600 for (i=0; i<B->rmap->n; i++) { 3601 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]); 3602 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); 3603 } 3604 } 3605 3606 B->preallocated = PETSC_TRUE; 3607 3608 b = (Mat_SeqAIJ*)B->data; 3609 3610 if (!skipallocation) { 3611 if (!b->imax) { 3612 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3613 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3614 } 3615 if (!b->ipre) { 3616 ierr = PetscMalloc1(B->rmap->n,&b->ipre);CHKERRQ(ierr); 3617 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3618 } 3619 if (!nnz) { 3620 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3621 else if (nz < 0) nz = 1; 3622 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3623 nz = nz*B->rmap->n; 3624 } else { 3625 nz = 0; 3626 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3627 } 3628 /* b->ilen will count nonzeros in each row so far. */ 3629 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3630 3631 /* allocate the matrix space */ 3632 /* FIXME: should B's old memory be unlogged? */ 3633 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3634 if (B->structure_only) { 3635 ierr = PetscMalloc1(nz,&b->j);CHKERRQ(ierr); 3636 ierr = PetscMalloc1(B->rmap->n+1,&b->i);CHKERRQ(ierr); 3637 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3638 } else { 3639 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3640 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3641 } 3642 b->i[0] = 0; 3643 for (i=1; i<B->rmap->n+1; i++) { 3644 b->i[i] = b->i[i-1] + b->imax[i-1]; 3645 } 3646 if (B->structure_only) { 3647 b->singlemalloc = PETSC_FALSE; 3648 b->free_a = PETSC_FALSE; 3649 } else { 3650 b->singlemalloc = PETSC_TRUE; 3651 b->free_a = PETSC_TRUE; 3652 } 3653 b->free_ij = PETSC_TRUE; 3654 } else { 3655 b->free_a = PETSC_FALSE; 3656 b->free_ij = PETSC_FALSE; 3657 } 3658 3659 if (b->ipre && nnz != b->ipre && b->imax) { 3660 /* reserve user-requested sparsity */ 3661 ierr = PetscMemcpy(b->ipre,b->imax,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3662 } 3663 3664 3665 b->nz = 0; 3666 b->maxnz = nz; 3667 B->info.nz_unneeded = (double)b->maxnz; 3668 if (realalloc) { 3669 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3670 } 3671 B->was_assembled = PETSC_FALSE; 3672 B->assembled = PETSC_FALSE; 3673 PetscFunctionReturn(0); 3674 } 3675 3676 3677 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 3678 { 3679 Mat_SeqAIJ *a; 3680 PetscInt i; 3681 PetscErrorCode ierr; 3682 3683 PetscFunctionBegin; 3684 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3685 a = (Mat_SeqAIJ*)A->data; 3686 /* if no saved info, we error out */ 3687 if (!a->ipre) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No saved preallocation info \n"); 3688 3689 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"); 3690 3691 ierr = PetscMemcpy(a->imax,a->ipre,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3692 ierr = PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3693 a->i[0] = 0; 3694 for (i=1; i<A->rmap->n+1; i++) { 3695 a->i[i] = a->i[i-1] + a->imax[i-1]; 3696 } 3697 A->preallocated = PETSC_TRUE; 3698 a->nz = 0; 3699 a->maxnz = a->i[A->rmap->n]; 3700 A->info.nz_unneeded = (double)a->maxnz; 3701 A->was_assembled = PETSC_FALSE; 3702 A->assembled = PETSC_FALSE; 3703 PetscFunctionReturn(0); 3704 } 3705 3706 /*@ 3707 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3708 3709 Input Parameters: 3710 + B - the matrix 3711 . i - the indices into j for the start of each row (starts with zero) 3712 . j - the column indices for each row (starts with zero) these must be sorted for each row 3713 - v - optional values in the matrix 3714 3715 Level: developer 3716 3717 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3718 3719 .keywords: matrix, aij, compressed row, sparse, sequential 3720 3721 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3722 @*/ 3723 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3724 { 3725 PetscErrorCode ierr; 3726 3727 PetscFunctionBegin; 3728 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3729 PetscValidType(B,1); 3730 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3731 PetscFunctionReturn(0); 3732 } 3733 3734 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3735 { 3736 PetscInt i; 3737 PetscInt m,n; 3738 PetscInt nz; 3739 PetscInt *nnz, nz_max = 0; 3740 PetscScalar *values; 3741 PetscErrorCode ierr; 3742 3743 PetscFunctionBegin; 3744 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3745 3746 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3747 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3748 3749 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3750 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3751 for (i = 0; i < m; i++) { 3752 nz = Ii[i+1]- Ii[i]; 3753 nz_max = PetscMax(nz_max, nz); 3754 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3755 nnz[i] = nz; 3756 } 3757 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3758 ierr = PetscFree(nnz);CHKERRQ(ierr); 3759 3760 if (v) { 3761 values = (PetscScalar*) v; 3762 } else { 3763 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3764 } 3765 3766 for (i = 0; i < m; i++) { 3767 nz = Ii[i+1] - Ii[i]; 3768 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3769 } 3770 3771 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3772 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3773 3774 if (!v) { 3775 ierr = PetscFree(values);CHKERRQ(ierr); 3776 } 3777 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3778 PetscFunctionReturn(0); 3779 } 3780 3781 #include <../src/mat/impls/dense/seq/dense.h> 3782 #include <petsc/private/kernels/petscaxpy.h> 3783 3784 /* 3785 Computes (B'*A')' since computing B*A directly is untenable 3786 3787 n p p 3788 ( ) ( ) ( ) 3789 m ( A ) * n ( B ) = m ( C ) 3790 ( ) ( ) ( ) 3791 3792 */ 3793 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3794 { 3795 PetscErrorCode ierr; 3796 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3797 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3798 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3799 PetscInt i,n,m,q,p; 3800 const PetscInt *ii,*idx; 3801 const PetscScalar *b,*a,*a_q; 3802 PetscScalar *c,*c_q; 3803 3804 PetscFunctionBegin; 3805 m = A->rmap->n; 3806 n = A->cmap->n; 3807 p = B->cmap->n; 3808 a = sub_a->v; 3809 b = sub_b->a; 3810 c = sub_c->v; 3811 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3812 3813 ii = sub_b->i; 3814 idx = sub_b->j; 3815 for (i=0; i<n; i++) { 3816 q = ii[i+1] - ii[i]; 3817 while (q-->0) { 3818 c_q = c + m*(*idx); 3819 a_q = a + m*i; 3820 PetscKernelAXPY(c_q,*b,a_q,m); 3821 idx++; 3822 b++; 3823 } 3824 } 3825 PetscFunctionReturn(0); 3826 } 3827 3828 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3829 { 3830 PetscErrorCode ierr; 3831 PetscInt m=A->rmap->n,n=B->cmap->n; 3832 Mat Cmat; 3833 3834 PetscFunctionBegin; 3835 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); 3836 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3837 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3838 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3839 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3840 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3841 3842 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3843 3844 *C = Cmat; 3845 PetscFunctionReturn(0); 3846 } 3847 3848 /* ----------------------------------------------------------------*/ 3849 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3850 { 3851 PetscErrorCode ierr; 3852 3853 PetscFunctionBegin; 3854 if (scall == MAT_INITIAL_MATRIX) { 3855 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3856 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3857 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3858 } 3859 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3860 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3861 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3862 PetscFunctionReturn(0); 3863 } 3864 3865 3866 /*MC 3867 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3868 based on compressed sparse row format. 3869 3870 Options Database Keys: 3871 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3872 3873 Level: beginner 3874 3875 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3876 M*/ 3877 3878 /*MC 3879 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3880 3881 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3882 and MATMPIAIJ otherwise. As a result, for single process communicators, 3883 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3884 for communicators controlling multiple processes. It is recommended that you call both of 3885 the above preallocation routines for simplicity. 3886 3887 Options Database Keys: 3888 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3889 3890 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 3891 enough exist. 3892 3893 Level: beginner 3894 3895 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3896 M*/ 3897 3898 /*MC 3899 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3900 3901 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3902 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3903 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3904 for communicators controlling multiple processes. It is recommended that you call both of 3905 the above preallocation routines for simplicity. 3906 3907 Options Database Keys: 3908 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3909 3910 Level: beginner 3911 3912 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3913 M*/ 3914 3915 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3916 #if defined(PETSC_HAVE_ELEMENTAL) 3917 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3918 #endif 3919 #if defined(PETSC_HAVE_HYPRE) 3920 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 3921 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 3922 #endif 3923 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3924 3925 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3926 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3927 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3928 #endif 3929 3930 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 3931 3932 /*@C 3933 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 3934 3935 Not Collective 3936 3937 Input Parameter: 3938 . mat - a MATSEQAIJ matrix 3939 3940 Output Parameter: 3941 . array - pointer to the data 3942 3943 Level: intermediate 3944 3945 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3946 @*/ 3947 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3948 { 3949 PetscErrorCode ierr; 3950 3951 PetscFunctionBegin; 3952 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3953 PetscFunctionReturn(0); 3954 } 3955 3956 /*@C 3957 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3958 3959 Not Collective 3960 3961 Input Parameter: 3962 . mat - a MATSEQAIJ matrix 3963 3964 Output Parameter: 3965 . nz - the maximum number of nonzeros in any row 3966 3967 Level: intermediate 3968 3969 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3970 @*/ 3971 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 3972 { 3973 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 3974 3975 PetscFunctionBegin; 3976 *nz = aij->rmax; 3977 PetscFunctionReturn(0); 3978 } 3979 3980 /*@C 3981 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 3982 3983 Not Collective 3984 3985 Input Parameters: 3986 . mat - a MATSEQAIJ matrix 3987 . array - pointer to the data 3988 3989 Level: intermediate 3990 3991 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3992 @*/ 3993 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3994 { 3995 PetscErrorCode ierr; 3996 3997 PetscFunctionBegin; 3998 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3999 PetscFunctionReturn(0); 4000 } 4001 4002 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4003 { 4004 Mat_SeqAIJ *b; 4005 PetscErrorCode ierr; 4006 PetscMPIInt size; 4007 4008 PetscFunctionBegin; 4009 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4010 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4011 4012 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4013 4014 B->data = (void*)b; 4015 4016 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4017 4018 b->row = 0; 4019 b->col = 0; 4020 b->icol = 0; 4021 b->reallocs = 0; 4022 b->ignorezeroentries = PETSC_FALSE; 4023 b->roworiented = PETSC_TRUE; 4024 b->nonew = 0; 4025 b->diag = 0; 4026 b->solve_work = 0; 4027 B->spptr = 0; 4028 b->saved_values = 0; 4029 b->idiag = 0; 4030 b->mdiag = 0; 4031 b->ssor_work = 0; 4032 b->omega = 1.0; 4033 b->fshift = 0.0; 4034 b->idiagvalid = PETSC_FALSE; 4035 b->ibdiagvalid = PETSC_FALSE; 4036 b->keepnonzeropattern = PETSC_FALSE; 4037 4038 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4039 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4040 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4041 4042 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4043 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4044 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4045 #endif 4046 4047 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4048 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4049 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4050 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4051 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4052 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4053 #if defined(PETSC_HAVE_MKL_SPARSE) 4054 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4055 #endif 4056 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4057 #if defined(PETSC_HAVE_ELEMENTAL) 4058 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4059 #endif 4060 #if defined(PETSC_HAVE_HYPRE) 4061 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 4062 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 4063 #endif 4064 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4065 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);CHKERRQ(ierr); 4066 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4067 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4068 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4069 ierr = PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);CHKERRQ(ierr); 4070 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4071 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4072 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4073 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4074 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4075 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4076 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4077 ierr = MatSeqAIJSetTypeFromOptions(B);CHKERRQ(ierr); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4078 PetscFunctionReturn(0); 4079 } 4080 4081 /* 4082 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4083 */ 4084 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4085 { 4086 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4087 PetscErrorCode ierr; 4088 PetscInt i,m = A->rmap->n; 4089 4090 PetscFunctionBegin; 4091 c = (Mat_SeqAIJ*)C->data; 4092 4093 C->factortype = A->factortype; 4094 c->row = 0; 4095 c->col = 0; 4096 c->icol = 0; 4097 c->reallocs = 0; 4098 4099 C->assembled = PETSC_TRUE; 4100 4101 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4102 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4103 4104 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4105 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4106 for (i=0; i<m; i++) { 4107 c->imax[i] = a->imax[i]; 4108 c->ilen[i] = a->ilen[i]; 4109 } 4110 4111 /* allocate the matrix space */ 4112 if (mallocmatspace) { 4113 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4114 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4115 4116 c->singlemalloc = PETSC_TRUE; 4117 4118 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4119 if (m > 0) { 4120 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4121 if (cpvalues == MAT_COPY_VALUES) { 4122 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4123 } else { 4124 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4125 } 4126 } 4127 } 4128 4129 c->ignorezeroentries = a->ignorezeroentries; 4130 c->roworiented = a->roworiented; 4131 c->nonew = a->nonew; 4132 if (a->diag) { 4133 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4134 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4135 for (i=0; i<m; i++) { 4136 c->diag[i] = a->diag[i]; 4137 } 4138 } else c->diag = 0; 4139 4140 c->solve_work = 0; 4141 c->saved_values = 0; 4142 c->idiag = 0; 4143 c->ssor_work = 0; 4144 c->keepnonzeropattern = a->keepnonzeropattern; 4145 c->free_a = PETSC_TRUE; 4146 c->free_ij = PETSC_TRUE; 4147 4148 c->rmax = a->rmax; 4149 c->nz = a->nz; 4150 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4151 C->preallocated = PETSC_TRUE; 4152 4153 c->compressedrow.use = a->compressedrow.use; 4154 c->compressedrow.nrows = a->compressedrow.nrows; 4155 if (a->compressedrow.use) { 4156 i = a->compressedrow.nrows; 4157 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4158 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4159 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4160 } else { 4161 c->compressedrow.use = PETSC_FALSE; 4162 c->compressedrow.i = NULL; 4163 c->compressedrow.rindex = NULL; 4164 } 4165 c->nonzerorowcnt = a->nonzerorowcnt; 4166 C->nonzerostate = A->nonzerostate; 4167 4168 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4169 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4170 PetscFunctionReturn(0); 4171 } 4172 4173 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4174 { 4175 PetscErrorCode ierr; 4176 4177 PetscFunctionBegin; 4178 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4179 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4180 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4181 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4182 } 4183 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4184 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4185 PetscFunctionReturn(0); 4186 } 4187 4188 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4189 { 4190 Mat_SeqAIJ *a; 4191 PetscErrorCode ierr; 4192 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4193 int fd; 4194 PetscMPIInt size; 4195 MPI_Comm comm; 4196 PetscInt bs = newMat->rmap->bs; 4197 4198 PetscFunctionBegin; 4199 /* force binary viewer to load .info file if it has not yet done so */ 4200 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4201 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4202 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4203 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4204 4205 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4206 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4207 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4208 if (bs < 0) bs = 1; 4209 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4210 4211 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4212 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4213 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4214 M = header[1]; N = header[2]; nz = header[3]; 4215 4216 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4217 4218 /* read in row lengths */ 4219 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4220 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4221 4222 /* check if sum of rowlengths is same as nz */ 4223 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4224 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); 4225 4226 /* set global size if not set already*/ 4227 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4228 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4229 } else { 4230 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4231 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4232 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4233 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4234 } 4235 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); 4236 } 4237 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4238 a = (Mat_SeqAIJ*)newMat->data; 4239 4240 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4241 4242 /* read in nonzero values */ 4243 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4244 4245 /* set matrix "i" values */ 4246 a->i[0] = 0; 4247 for (i=1; i<= M; i++) { 4248 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4249 a->ilen[i-1] = rowlengths[i-1]; 4250 } 4251 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4252 4253 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4254 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4255 PetscFunctionReturn(0); 4256 } 4257 4258 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4259 { 4260 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4261 PetscErrorCode ierr; 4262 #if defined(PETSC_USE_COMPLEX) 4263 PetscInt k; 4264 #endif 4265 4266 PetscFunctionBegin; 4267 /* If the matrix dimensions are not equal,or no of nonzeros */ 4268 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4269 *flg = PETSC_FALSE; 4270 PetscFunctionReturn(0); 4271 } 4272 4273 /* if the a->i are the same */ 4274 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4275 if (!*flg) PetscFunctionReturn(0); 4276 4277 /* if a->j are the same */ 4278 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4279 if (!*flg) PetscFunctionReturn(0); 4280 4281 /* if a->a are the same */ 4282 #if defined(PETSC_USE_COMPLEX) 4283 for (k=0; k<a->nz; k++) { 4284 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4285 *flg = PETSC_FALSE; 4286 PetscFunctionReturn(0); 4287 } 4288 } 4289 #else 4290 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4291 #endif 4292 PetscFunctionReturn(0); 4293 } 4294 4295 /*@ 4296 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4297 provided by the user. 4298 4299 Collective on MPI_Comm 4300 4301 Input Parameters: 4302 + comm - must be an MPI communicator of size 1 4303 . m - number of rows 4304 . n - number of columns 4305 . i - row indices 4306 . j - column indices 4307 - a - matrix values 4308 4309 Output Parameter: 4310 . mat - the matrix 4311 4312 Level: intermediate 4313 4314 Notes: 4315 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4316 once the matrix is destroyed and not before 4317 4318 You cannot set new nonzero locations into this matrix, that will generate an error. 4319 4320 The i and j indices are 0 based 4321 4322 The format which is used for the sparse matrix input, is equivalent to a 4323 row-major ordering.. i.e for the following matrix, the input data expected is 4324 as shown 4325 4326 $ 1 0 0 4327 $ 2 0 3 4328 $ 4 5 6 4329 $ 4330 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4331 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4332 $ v = {1,2,3,4,5,6} [size = 6] 4333 4334 4335 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4336 4337 @*/ 4338 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4339 { 4340 PetscErrorCode ierr; 4341 PetscInt ii; 4342 Mat_SeqAIJ *aij; 4343 #if defined(PETSC_USE_DEBUG) 4344 PetscInt jj; 4345 #endif 4346 4347 PetscFunctionBegin; 4348 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4349 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4350 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4351 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4352 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4353 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4354 aij = (Mat_SeqAIJ*)(*mat)->data; 4355 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4356 4357 aij->i = i; 4358 aij->j = j; 4359 aij->a = a; 4360 aij->singlemalloc = PETSC_FALSE; 4361 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4362 aij->free_a = PETSC_FALSE; 4363 aij->free_ij = PETSC_FALSE; 4364 4365 for (ii=0; ii<m; ii++) { 4366 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4367 #if defined(PETSC_USE_DEBUG) 4368 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]); 4369 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4370 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); 4371 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); 4372 } 4373 #endif 4374 } 4375 #if defined(PETSC_USE_DEBUG) 4376 for (ii=0; ii<aij->i[m]; ii++) { 4377 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4378 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]); 4379 } 4380 #endif 4381 4382 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4383 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4384 PetscFunctionReturn(0); 4385 } 4386 /*@C 4387 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4388 provided by the user. 4389 4390 Collective on MPI_Comm 4391 4392 Input Parameters: 4393 + comm - must be an MPI communicator of size 1 4394 . m - number of rows 4395 . n - number of columns 4396 . i - row indices 4397 . j - column indices 4398 . a - matrix values 4399 . nz - number of nonzeros 4400 - idx - 0 or 1 based 4401 4402 Output Parameter: 4403 . mat - the matrix 4404 4405 Level: intermediate 4406 4407 Notes: 4408 The i and j indices are 0 based 4409 4410 The format which is used for the sparse matrix input, is equivalent to a 4411 row-major ordering.. i.e for the following matrix, the input data expected is 4412 as shown: 4413 4414 1 0 0 4415 2 0 3 4416 4 5 6 4417 4418 i = {0,1,1,2,2,2} 4419 j = {0,0,2,0,1,2} 4420 v = {1,2,3,4,5,6} 4421 4422 4423 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4424 4425 @*/ 4426 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4427 { 4428 PetscErrorCode ierr; 4429 PetscInt ii, *nnz, one = 1,row,col; 4430 4431 4432 PetscFunctionBegin; 4433 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4434 for (ii = 0; ii < nz; ii++) { 4435 nnz[i[ii] - !!idx] += 1; 4436 } 4437 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4438 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4439 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4440 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4441 for (ii = 0; ii < nz; ii++) { 4442 if (idx) { 4443 row = i[ii] - 1; 4444 col = j[ii] - 1; 4445 } else { 4446 row = i[ii]; 4447 col = j[ii]; 4448 } 4449 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4450 } 4451 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4452 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4453 ierr = PetscFree(nnz);CHKERRQ(ierr); 4454 PetscFunctionReturn(0); 4455 } 4456 4457 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4458 { 4459 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4460 PetscErrorCode ierr; 4461 4462 PetscFunctionBegin; 4463 a->idiagvalid = PETSC_FALSE; 4464 a->ibdiagvalid = PETSC_FALSE; 4465 4466 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4467 PetscFunctionReturn(0); 4468 } 4469 4470 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4471 { 4472 PetscErrorCode ierr; 4473 PetscMPIInt size; 4474 4475 PetscFunctionBegin; 4476 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4477 if (size == 1) { 4478 if (scall == MAT_INITIAL_MATRIX) { 4479 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 4480 } else { 4481 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4482 } 4483 } else { 4484 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4485 } 4486 PetscFunctionReturn(0); 4487 } 4488 4489 /* 4490 Permute A into C's *local* index space using rowemb,colemb. 4491 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4492 of [0,m), colemb is in [0,n). 4493 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4494 */ 4495 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4496 { 4497 /* If making this function public, change the error returned in this function away from _PLIB. */ 4498 PetscErrorCode ierr; 4499 Mat_SeqAIJ *Baij; 4500 PetscBool seqaij; 4501 PetscInt m,n,*nz,i,j,count; 4502 PetscScalar v; 4503 const PetscInt *rowindices,*colindices; 4504 4505 PetscFunctionBegin; 4506 if (!B) PetscFunctionReturn(0); 4507 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4508 ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4509 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4510 if (rowemb) { 4511 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4512 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); 4513 } else { 4514 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4515 } 4516 if (colemb) { 4517 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4518 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); 4519 } else { 4520 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4521 } 4522 4523 Baij = (Mat_SeqAIJ*)(B->data); 4524 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4525 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4526 for (i=0; i<B->rmap->n; i++) { 4527 nz[i] = Baij->i[i+1] - Baij->i[i]; 4528 } 4529 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4530 ierr = PetscFree(nz);CHKERRQ(ierr); 4531 } 4532 if (pattern == SUBSET_NONZERO_PATTERN) { 4533 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4534 } 4535 count = 0; 4536 rowindices = NULL; 4537 colindices = NULL; 4538 if (rowemb) { 4539 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4540 } 4541 if (colemb) { 4542 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4543 } 4544 for (i=0; i<B->rmap->n; i++) { 4545 PetscInt row; 4546 row = i; 4547 if (rowindices) row = rowindices[i]; 4548 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4549 PetscInt col; 4550 col = Baij->j[count]; 4551 if (colindices) col = colindices[col]; 4552 v = Baij->a[count]; 4553 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4554 ++count; 4555 } 4556 } 4557 /* FIXME: set C's nonzerostate correctly. */ 4558 /* Assembly for C is necessary. */ 4559 C->preallocated = PETSC_TRUE; 4560 C->assembled = PETSC_TRUE; 4561 C->was_assembled = PETSC_FALSE; 4562 PetscFunctionReturn(0); 4563 } 4564 4565 PetscFunctionList MatSeqAIJList = NULL; 4566 4567 /*@C 4568 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 4569 4570 Collective on Mat 4571 4572 Input Parameters: 4573 + mat - the matrix object 4574 - matype - matrix type 4575 4576 Options Database Key: 4577 . -mat_seqai_type <method> - for example seqaijcrl 4578 4579 4580 Level: intermediate 4581 4582 .keywords: Mat, MatType, set, method 4583 4584 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 4585 @*/ 4586 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 4587 { 4588 PetscErrorCode ierr,(*r)(Mat,const MatType,MatReuse,Mat*); 4589 PetscBool sametype; 4590 4591 PetscFunctionBegin; 4592 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 4593 ierr = PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr); 4594 if (sametype) PetscFunctionReturn(0); 4595 4596 ierr = PetscFunctionListFind(MatSeqAIJList,matype,&r);CHKERRQ(ierr); 4597 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 4598 ierr = (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);CHKERRQ(ierr); 4599 PetscFunctionReturn(0); 4600 } 4601 4602 4603 /*@C 4604 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 4605 4606 Not Collective 4607 4608 Input Parameters: 4609 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 4610 - function - routine to convert to subtype 4611 4612 Notes: 4613 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 4614 4615 4616 Then, your matrix can be chosen with the procedural interface at runtime via the option 4617 $ -mat_seqaij_type my_mat 4618 4619 Level: advanced 4620 4621 .keywords: Mat, register 4622 4623 .seealso: MatSeqAIJRegisterAll() 4624 4625 4626 Level: advanced 4627 @*/ 4628 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 4629 { 4630 PetscErrorCode ierr; 4631 4632 PetscFunctionBegin; 4633 ierr = PetscFunctionListAdd(&MatSeqAIJList,sname,function);CHKERRQ(ierr); 4634 PetscFunctionReturn(0); 4635 } 4636 4637 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 4638 4639 /*@C 4640 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 4641 4642 Not Collective 4643 4644 Level: advanced 4645 4646 Developers Note: CUSP and CUSPARSE do not yet support the MatConvert_SeqAIJ..() paradigm and thus cannot be registered here 4647 4648 .keywords: KSP, register, all 4649 4650 .seealso: MatRegisterAll(), MatSeqAIJRegister() 4651 @*/ 4652 PetscErrorCode MatSeqAIJRegisterAll(void) 4653 { 4654 PetscErrorCode ierr; 4655 4656 PetscFunctionBegin; 4657 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 4658 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 4659 4660 ierr = MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4661 ierr = MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4662 #if defined(PETSC_HAVE_MKL_SPARSE) 4663 ierr = MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4664 #endif 4665 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 4666 ierr = MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 4667 #endif 4668 PetscFunctionReturn(0); 4669 } 4670 4671 /* 4672 Special version for direct calls from Fortran 4673 */ 4674 #include <petsc/private/fortranimpl.h> 4675 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4676 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4677 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4678 #define matsetvaluesseqaij_ matsetvaluesseqaij 4679 #endif 4680 4681 /* Change these macros so can be used in void function */ 4682 #undef CHKERRQ 4683 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4684 #undef SETERRQ2 4685 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4686 #undef SETERRQ3 4687 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4688 4689 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) 4690 { 4691 Mat A = *AA; 4692 PetscInt m = *mm, n = *nn; 4693 InsertMode is = *isis; 4694 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4695 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4696 PetscInt *imax,*ai,*ailen; 4697 PetscErrorCode ierr; 4698 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4699 MatScalar *ap,value,*aa; 4700 PetscBool ignorezeroentries = a->ignorezeroentries; 4701 PetscBool roworiented = a->roworiented; 4702 4703 PetscFunctionBegin; 4704 MatCheckPreallocated(A,1); 4705 imax = a->imax; 4706 ai = a->i; 4707 ailen = a->ilen; 4708 aj = a->j; 4709 aa = a->a; 4710 4711 for (k=0; k<m; k++) { /* loop over added rows */ 4712 row = im[k]; 4713 if (row < 0) continue; 4714 #if defined(PETSC_USE_DEBUG) 4715 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4716 #endif 4717 rp = aj + ai[row]; ap = aa + ai[row]; 4718 rmax = imax[row]; nrow = ailen[row]; 4719 low = 0; 4720 high = nrow; 4721 for (l=0; l<n; l++) { /* loop over added columns */ 4722 if (in[l] < 0) continue; 4723 #if defined(PETSC_USE_DEBUG) 4724 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4725 #endif 4726 col = in[l]; 4727 if (roworiented) value = v[l + k*n]; 4728 else value = v[k + l*m]; 4729 4730 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4731 4732 if (col <= lastcol) low = 0; 4733 else high = nrow; 4734 lastcol = col; 4735 while (high-low > 5) { 4736 t = (low+high)/2; 4737 if (rp[t] > col) high = t; 4738 else low = t; 4739 } 4740 for (i=low; i<high; i++) { 4741 if (rp[i] > col) break; 4742 if (rp[i] == col) { 4743 if (is == ADD_VALUES) ap[i] += value; 4744 else ap[i] = value; 4745 goto noinsert; 4746 } 4747 } 4748 if (value == 0.0 && ignorezeroentries) goto noinsert; 4749 if (nonew == 1) goto noinsert; 4750 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4751 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4752 N = nrow++ - 1; a->nz++; high++; 4753 /* shift up all the later entries in this row */ 4754 for (ii=N; ii>=i; ii--) { 4755 rp[ii+1] = rp[ii]; 4756 ap[ii+1] = ap[ii]; 4757 } 4758 rp[i] = col; 4759 ap[i] = value; 4760 A->nonzerostate++; 4761 noinsert:; 4762 low = i + 1; 4763 } 4764 ailen[row] = nrow; 4765 } 4766 PetscFunctionReturnVoid(); 4767 } 4768 4769