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