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