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