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