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