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