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