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