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