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