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