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