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