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