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