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