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