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