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