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 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2045 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2046 ierr = MatSetBlockSizes(C,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 2047 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2048 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2049 } 2050 c = (Mat_SeqAIJ*)C->data; 2051 2052 /* loop over rows inserting into submatrix */ 2053 a_new = c->a; 2054 j_new = c->j; 2055 i_new = c->i; 2056 2057 for (i=0; i<nrows; i++) { 2058 ii = starts[i]; 2059 lensi = lens[i]; 2060 for (k=0; k<lensi; k++) { 2061 *j_new++ = aj[ii+k] - first; 2062 } 2063 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2064 a_new += lensi; 2065 i_new[i+1] = i_new[i] + lensi; 2066 c->ilen[i] = lensi; 2067 } 2068 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2069 } else { 2070 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2071 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 2072 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 2073 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2074 for (i=0; i<ncols; i++) { 2075 #if defined(PETSC_USE_DEBUG) 2076 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); 2077 #endif 2078 smap[icol[i]] = i+1; 2079 } 2080 2081 /* determine lens of each row */ 2082 for (i=0; i<nrows; i++) { 2083 kstart = ai[irow[i]]; 2084 kend = kstart + a->ilen[irow[i]]; 2085 lens[i] = 0; 2086 for (k=kstart; k<kend; k++) { 2087 if (smap[aj[k]]) { 2088 lens[i]++; 2089 } 2090 } 2091 } 2092 /* Create and fill new matrix */ 2093 if (scall == MAT_REUSE_MATRIX) { 2094 PetscBool equal; 2095 2096 c = (Mat_SeqAIJ *)((*B)->data); 2097 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2098 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2099 if (!equal) { 2100 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2101 } 2102 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2103 C = *B; 2104 } else { 2105 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2106 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2107 ierr = MatSetBlockSizes(C,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 2108 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2109 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2110 } 2111 c = (Mat_SeqAIJ *)(C->data); 2112 for (i=0; i<nrows; i++) { 2113 row = irow[i]; 2114 kstart = ai[row]; 2115 kend = kstart + a->ilen[row]; 2116 mat_i = c->i[i]; 2117 mat_j = c->j + mat_i; 2118 mat_a = c->a + mat_i; 2119 mat_ilen = c->ilen + i; 2120 for (k=kstart; k<kend; k++) { 2121 if ((tcol=smap[a->j[k]])) { 2122 *mat_j++ = tcol - 1; 2123 *mat_a++ = a->a[k]; 2124 (*mat_ilen)++; 2125 2126 } 2127 } 2128 } 2129 /* Free work space */ 2130 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2131 ierr = PetscFree(smap);CHKERRQ(ierr); 2132 ierr = PetscFree(lens);CHKERRQ(ierr); 2133 } 2134 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2135 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2136 2137 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2138 *B = C; 2139 PetscFunctionReturn(0); 2140 } 2141 2142 #undef __FUNCT__ 2143 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2144 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat* subMat) 2145 { 2146 PetscErrorCode ierr; 2147 Mat B; 2148 2149 PetscFunctionBegin; 2150 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2151 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2152 ierr = MatSetBlockSizes(B,mat->rmap->bs,mat->cmap->bs); CHKERRQ(ierr); 2153 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2154 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2155 *subMat = B; 2156 PetscFunctionReturn(0); 2157 } 2158 2159 #undef __FUNCT__ 2160 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2161 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2162 { 2163 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2164 PetscErrorCode ierr; 2165 Mat outA; 2166 PetscBool row_identity,col_identity; 2167 2168 PetscFunctionBegin; 2169 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2170 2171 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2172 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2173 2174 outA = inA; 2175 outA->factortype = MAT_FACTOR_LU; 2176 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2177 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2178 a->row = row; 2179 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2180 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2181 a->col = col; 2182 2183 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2184 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2185 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2186 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 2187 2188 if (!a->solve_work) { /* this matrix may have been factored before */ 2189 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 2190 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2191 } 2192 2193 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2194 if (row_identity && col_identity) { 2195 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2196 } else { 2197 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2198 } 2199 PetscFunctionReturn(0); 2200 } 2201 2202 #undef __FUNCT__ 2203 #define __FUNCT__ "MatScale_SeqAIJ" 2204 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2205 { 2206 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2207 PetscScalar oalpha = alpha; 2208 PetscErrorCode ierr; 2209 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 2210 2211 PetscFunctionBegin; 2212 BLASscal_(&bnz,&oalpha,a->a,&one); 2213 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2214 a->idiagvalid = PETSC_FALSE; 2215 a->ibdiagvalid = PETSC_FALSE; 2216 PetscFunctionReturn(0); 2217 } 2218 2219 #undef __FUNCT__ 2220 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2221 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2222 { 2223 PetscErrorCode ierr; 2224 PetscInt i; 2225 2226 PetscFunctionBegin; 2227 if (scall == MAT_INITIAL_MATRIX) { 2228 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 2229 } 2230 2231 for (i=0; i<n; i++) { 2232 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2233 } 2234 PetscFunctionReturn(0); 2235 } 2236 2237 #undef __FUNCT__ 2238 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2239 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2240 { 2241 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2242 PetscErrorCode ierr; 2243 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2244 const PetscInt *idx; 2245 PetscInt start,end,*ai,*aj; 2246 PetscBT table; 2247 2248 PetscFunctionBegin; 2249 m = A->rmap->n; 2250 ai = a->i; 2251 aj = a->j; 2252 2253 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2254 2255 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 2256 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2257 2258 for (i=0; i<is_max; i++) { 2259 /* Initialize the two local arrays */ 2260 isz = 0; 2261 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2262 2263 /* Extract the indices, assume there can be duplicate entries */ 2264 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2265 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2266 2267 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2268 for (j=0; j<n ; ++j){ 2269 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 2270 } 2271 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2272 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2273 2274 k = 0; 2275 for (j=0; j<ov; j++){ /* for each overlap */ 2276 n = isz; 2277 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 2278 row = nidx[k]; 2279 start = ai[row]; 2280 end = ai[row+1]; 2281 for (l = start; l<end ; l++){ 2282 val = aj[l] ; 2283 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2284 } 2285 } 2286 } 2287 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2288 } 2289 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2290 ierr = PetscFree(nidx);CHKERRQ(ierr); 2291 PetscFunctionReturn(0); 2292 } 2293 2294 /* -------------------------------------------------------------- */ 2295 #undef __FUNCT__ 2296 #define __FUNCT__ "MatPermute_SeqAIJ" 2297 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2298 { 2299 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2300 PetscErrorCode ierr; 2301 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2302 const PetscInt *row,*col; 2303 PetscInt *cnew,j,*lens; 2304 IS icolp,irowp; 2305 PetscInt *cwork = PETSC_NULL; 2306 PetscScalar *vwork = PETSC_NULL; 2307 2308 PetscFunctionBegin; 2309 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2310 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2311 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2312 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2313 2314 /* determine lengths of permuted rows */ 2315 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2316 for (i=0; i<m; i++) { 2317 lens[row[i]] = a->i[i+1] - a->i[i]; 2318 } 2319 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2320 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2321 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 2322 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2323 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2324 ierr = PetscFree(lens);CHKERRQ(ierr); 2325 2326 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2327 for (i=0; i<m; i++) { 2328 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2329 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2330 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2331 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2332 } 2333 ierr = PetscFree(cnew);CHKERRQ(ierr); 2334 (*B)->assembled = PETSC_FALSE; 2335 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2336 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2337 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2338 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2339 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2340 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2341 PetscFunctionReturn(0); 2342 } 2343 2344 #undef __FUNCT__ 2345 #define __FUNCT__ "MatCopy_SeqAIJ" 2346 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2347 { 2348 PetscErrorCode ierr; 2349 2350 PetscFunctionBegin; 2351 /* If the two matrices have the same copy implementation, use fast copy. */ 2352 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2353 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2354 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2355 2356 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"); 2357 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2358 } else { 2359 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2360 } 2361 PetscFunctionReturn(0); 2362 } 2363 2364 #undef __FUNCT__ 2365 #define __FUNCT__ "MatSetUp_SeqAIJ" 2366 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2367 { 2368 PetscErrorCode ierr; 2369 2370 PetscFunctionBegin; 2371 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2372 PetscFunctionReturn(0); 2373 } 2374 2375 #undef __FUNCT__ 2376 #define __FUNCT__ "MatGetArray_SeqAIJ" 2377 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2378 { 2379 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2380 PetscFunctionBegin; 2381 *array = a->a; 2382 PetscFunctionReturn(0); 2383 } 2384 2385 #undef __FUNCT__ 2386 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2387 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2388 { 2389 PetscFunctionBegin; 2390 PetscFunctionReturn(0); 2391 } 2392 2393 #undef __FUNCT__ 2394 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2395 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2396 { 2397 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2398 PetscErrorCode ierr; 2399 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2400 PetscScalar dx,*y,*xx,*w3_array; 2401 PetscScalar *vscale_array; 2402 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2403 Vec w1,w2,w3; 2404 void *fctx = coloring->fctx; 2405 PetscBool flg = PETSC_FALSE; 2406 2407 PetscFunctionBegin; 2408 if (!coloring->w1) { 2409 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2410 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2411 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2412 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2413 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2414 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2415 } 2416 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2417 2418 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2419 ierr = PetscOptionsGetBool(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2420 if (flg) { 2421 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2422 } else { 2423 PetscBool assembled; 2424 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2425 if (assembled) { 2426 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2427 } 2428 } 2429 2430 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2431 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2432 2433 /* 2434 This is a horrible, horrible, hack. 2435 */ 2436 if (coloring->F) { 2437 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2438 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2439 if (m1 != m2) { 2440 coloring->F = 0; 2441 } 2442 } 2443 2444 if (coloring->F) { 2445 w1 = coloring->F; 2446 coloring->F = 0; 2447 } else { 2448 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2449 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2450 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2451 } 2452 2453 /* 2454 Compute all the scale factors and share with other processors 2455 */ 2456 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2457 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2458 for (k=0; k<coloring->ncolors; k++) { 2459 /* 2460 Loop over each column associated with color adding the 2461 perturbation to the vector w3. 2462 */ 2463 for (l=0; l<coloring->ncolumns[k]; l++) { 2464 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2465 dx = xx[col]; 2466 if (dx == 0.0) dx = 1.0; 2467 #if !defined(PETSC_USE_COMPLEX) 2468 if (dx < umin && dx >= 0.0) dx = umin; 2469 else if (dx < 0.0 && dx > -umin) dx = -umin; 2470 #else 2471 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2472 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2473 #endif 2474 dx *= epsilon; 2475 vscale_array[col] = 1.0/dx; 2476 } 2477 } 2478 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2479 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2480 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2481 2482 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2483 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2484 2485 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2486 else vscaleforrow = coloring->columnsforrow; 2487 2488 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2489 /* 2490 Loop over each color 2491 */ 2492 for (k=0; k<coloring->ncolors; k++) { 2493 coloring->currentcolor = k; 2494 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2495 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2496 /* 2497 Loop over each column associated with color adding the 2498 perturbation to the vector w3. 2499 */ 2500 for (l=0; l<coloring->ncolumns[k]; l++) { 2501 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2502 dx = xx[col]; 2503 if (dx == 0.0) dx = 1.0; 2504 #if !defined(PETSC_USE_COMPLEX) 2505 if (dx < umin && dx >= 0.0) dx = umin; 2506 else if (dx < 0.0 && dx > -umin) dx = -umin; 2507 #else 2508 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2509 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2510 #endif 2511 dx *= epsilon; 2512 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2513 w3_array[col] += dx; 2514 } 2515 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2516 2517 /* 2518 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2519 */ 2520 2521 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2522 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2523 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2524 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2525 2526 /* 2527 Loop over rows of vector, putting results into Jacobian matrix 2528 */ 2529 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2530 for (l=0; l<coloring->nrows[k]; l++) { 2531 row = coloring->rows[k][l]; 2532 col = coloring->columnsforrow[k][l]; 2533 y[row] *= vscale_array[vscaleforrow[k][l]]; 2534 srow = row + start; 2535 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2536 } 2537 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2538 } 2539 coloring->currentcolor = k; 2540 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2541 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2542 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2543 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2544 PetscFunctionReturn(0); 2545 } 2546 2547 /* 2548 Computes the number of nonzeros per row needed for preallocation when X and Y 2549 have different nonzero structure. 2550 */ 2551 #undef __FUNCT__ 2552 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2553 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt* nnz) 2554 { 2555 PetscInt i,m=Y->rmap->N; 2556 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2557 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2558 const PetscInt *xi = x->i,*yi = y->i; 2559 2560 PetscFunctionBegin; 2561 /* Set the number of nonzeros in the new matrix */ 2562 for(i=0; i<m; i++) { 2563 PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i]; 2564 const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i]; 2565 nnz[i] = 0; 2566 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2567 for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */ 2568 if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */ 2569 nnz[i]++; 2570 } 2571 for (; k<nzy; k++) nnz[i]++; 2572 } 2573 PetscFunctionReturn(0); 2574 } 2575 2576 #undef __FUNCT__ 2577 #define __FUNCT__ "MatAXPY_SeqAIJ" 2578 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2579 { 2580 PetscErrorCode ierr; 2581 PetscInt i; 2582 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2583 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2584 2585 PetscFunctionBegin; 2586 if (str == SAME_NONZERO_PATTERN) { 2587 PetscScalar alpha = a; 2588 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2589 y->idiagvalid = PETSC_FALSE; 2590 y->ibdiagvalid = PETSC_FALSE; 2591 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2592 if (y->xtoy && y->XtoY != X) { 2593 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2594 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 2595 } 2596 if (!y->xtoy) { /* get xtoy */ 2597 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2598 y->XtoY = X; 2599 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2600 } 2601 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2602 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2603 } else { 2604 Mat B; 2605 PetscInt *nnz; 2606 ierr = PetscMalloc(Y->rmap->N*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 2607 ierr = MatCreate(((PetscObject)Y)->comm,&B);CHKERRQ(ierr); 2608 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2609 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2610 ierr = MatSetBlockSizes(B,Y->rmap->bs,Y->cmap->bs);CHKERRQ(ierr); 2611 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2612 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2613 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2614 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2615 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2616 ierr = PetscFree(nnz);CHKERRQ(ierr); 2617 } 2618 PetscFunctionReturn(0); 2619 } 2620 2621 #undef __FUNCT__ 2622 #define __FUNCT__ "MatConjugate_SeqAIJ" 2623 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2624 { 2625 #if defined(PETSC_USE_COMPLEX) 2626 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2627 PetscInt i,nz; 2628 PetscScalar *a; 2629 2630 PetscFunctionBegin; 2631 nz = aij->nz; 2632 a = aij->a; 2633 for (i=0; i<nz; i++) { 2634 a[i] = PetscConj(a[i]); 2635 } 2636 #else 2637 PetscFunctionBegin; 2638 #endif 2639 PetscFunctionReturn(0); 2640 } 2641 2642 #undef __FUNCT__ 2643 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2644 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2645 { 2646 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2647 PetscErrorCode ierr; 2648 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2649 PetscReal atmp; 2650 PetscScalar *x; 2651 MatScalar *aa; 2652 2653 PetscFunctionBegin; 2654 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2655 aa = a->a; 2656 ai = a->i; 2657 aj = a->j; 2658 2659 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2660 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2661 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2662 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2663 for (i=0; i<m; i++) { 2664 ncols = ai[1] - ai[0]; ai++; 2665 x[i] = 0.0; 2666 for (j=0; j<ncols; j++){ 2667 atmp = PetscAbsScalar(*aa); 2668 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2669 aa++; aj++; 2670 } 2671 } 2672 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2673 PetscFunctionReturn(0); 2674 } 2675 2676 #undef __FUNCT__ 2677 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2678 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2679 { 2680 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2681 PetscErrorCode ierr; 2682 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2683 PetscScalar *x; 2684 MatScalar *aa; 2685 2686 PetscFunctionBegin; 2687 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2688 aa = a->a; 2689 ai = a->i; 2690 aj = a->j; 2691 2692 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2693 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2694 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2695 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2696 for (i=0; i<m; i++) { 2697 ncols = ai[1] - ai[0]; ai++; 2698 if (ncols == A->cmap->n) { /* row is dense */ 2699 x[i] = *aa; if (idx) idx[i] = 0; 2700 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2701 x[i] = 0.0; 2702 if (idx) { 2703 idx[i] = 0; /* in case ncols is zero */ 2704 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2705 if (aj[j] > j) { 2706 idx[i] = j; 2707 break; 2708 } 2709 } 2710 } 2711 } 2712 for (j=0; j<ncols; j++){ 2713 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2714 aa++; aj++; 2715 } 2716 } 2717 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2718 PetscFunctionReturn(0); 2719 } 2720 2721 #undef __FUNCT__ 2722 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2723 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2724 { 2725 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2726 PetscErrorCode ierr; 2727 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2728 PetscReal atmp; 2729 PetscScalar *x; 2730 MatScalar *aa; 2731 2732 PetscFunctionBegin; 2733 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2734 aa = a->a; 2735 ai = a->i; 2736 aj = a->j; 2737 2738 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2739 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2740 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2741 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); 2742 for (i=0; i<m; i++) { 2743 ncols = ai[1] - ai[0]; ai++; 2744 if (ncols) { 2745 /* Get first nonzero */ 2746 for(j = 0; j < ncols; j++) { 2747 atmp = PetscAbsScalar(aa[j]); 2748 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2749 } 2750 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2751 } else { 2752 x[i] = 0.0; if (idx) idx[i] = 0; 2753 } 2754 for(j = 0; j < ncols; j++) { 2755 atmp = PetscAbsScalar(*aa); 2756 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2757 aa++; aj++; 2758 } 2759 } 2760 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2761 PetscFunctionReturn(0); 2762 } 2763 2764 #undef __FUNCT__ 2765 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2766 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2767 { 2768 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2769 PetscErrorCode ierr; 2770 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2771 PetscScalar *x; 2772 MatScalar *aa; 2773 2774 PetscFunctionBegin; 2775 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2776 aa = a->a; 2777 ai = a->i; 2778 aj = a->j; 2779 2780 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2781 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2782 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2783 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2784 for (i=0; i<m; i++) { 2785 ncols = ai[1] - ai[0]; ai++; 2786 if (ncols == A->cmap->n) { /* row is dense */ 2787 x[i] = *aa; if (idx) idx[i] = 0; 2788 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2789 x[i] = 0.0; 2790 if (idx) { /* find first implicit 0.0 in the row */ 2791 idx[i] = 0; /* in case ncols is zero */ 2792 for (j=0;j<ncols;j++) { 2793 if (aj[j] > j) { 2794 idx[i] = j; 2795 break; 2796 } 2797 } 2798 } 2799 } 2800 for (j=0; j<ncols; j++){ 2801 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2802 aa++; aj++; 2803 } 2804 } 2805 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2806 PetscFunctionReturn(0); 2807 } 2808 2809 #include <petscblaslapack.h> 2810 #include <../src/mat/blockinvert.h> 2811 2812 #undef __FUNCT__ 2813 #define __FUNCT__ "MatInvertBlockDiagonal_SeqAIJ" 2814 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,PetscScalar **values) 2815 { 2816 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2817 PetscErrorCode ierr; 2818 PetscInt i,bs = A->rmap->bs,mbs = A->rmap->n/A->rmap->bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2819 MatScalar *diag,work[25],*v_work; 2820 PetscReal shift = 0.0; 2821 2822 PetscFunctionBegin; 2823 if (a->ibdiagvalid) { 2824 if (values) *values = a->ibdiag; 2825 PetscFunctionReturn(0); 2826 } 2827 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2828 if (!a->ibdiag) { 2829 ierr = PetscMalloc(bs2*mbs*sizeof(PetscScalar),&a->ibdiag);CHKERRQ(ierr); 2830 ierr = PetscLogObjectMemory(A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2831 } 2832 diag = a->ibdiag; 2833 if (values) *values = a->ibdiag; 2834 /* factor and invert each block */ 2835 switch (bs){ 2836 case 1: 2837 for (i=0; i<mbs; i++) { 2838 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 2839 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 2840 } 2841 break; 2842 case 2: 2843 for (i=0; i<mbs; i++) { 2844 ij[0] = 2*i; ij[1] = 2*i + 1; 2845 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 2846 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift);CHKERRQ(ierr); 2847 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 2848 diag += 4; 2849 } 2850 break; 2851 case 3: 2852 for (i=0; i<mbs; i++) { 2853 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 2854 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 2855 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 2856 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 2857 diag += 9; 2858 } 2859 break; 2860 case 4: 2861 for (i=0; i<mbs; i++) { 2862 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 2863 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 2864 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 2865 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 2866 diag += 16; 2867 } 2868 break; 2869 case 5: 2870 for (i=0; i<mbs; i++) { 2871 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 2872 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 2873 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 2874 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 2875 diag += 25; 2876 } 2877 break; 2878 case 6: 2879 for (i=0; i<mbs; i++) { 2880 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; 2881 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 2882 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 2883 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 2884 diag += 36; 2885 } 2886 break; 2887 case 7: 2888 for (i=0; i<mbs; i++) { 2889 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; 2890 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 2891 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 2892 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 2893 diag += 49; 2894 } 2895 break; 2896 default: 2897 ierr = PetscMalloc3(bs,MatScalar,&v_work,bs,PetscInt,&v_pivots,bs,PetscInt,&IJ);CHKERRQ(ierr); 2898 for (i=0; i<mbs; i++) { 2899 for (j=0; j<bs; j++) { 2900 IJ[j] = bs*i + j; 2901 } 2902 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 2903 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 2904 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 2905 diag += bs2; 2906 } 2907 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 2908 } 2909 a->ibdiagvalid = PETSC_TRUE; 2910 PetscFunctionReturn(0); 2911 } 2912 2913 extern PetscErrorCode MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2914 /* -------------------------------------------------------------------*/ 2915 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2916 MatGetRow_SeqAIJ, 2917 MatRestoreRow_SeqAIJ, 2918 MatMult_SeqAIJ, 2919 /* 4*/ MatMultAdd_SeqAIJ, 2920 MatMultTranspose_SeqAIJ, 2921 MatMultTransposeAdd_SeqAIJ, 2922 0, 2923 0, 2924 0, 2925 /*10*/ 0, 2926 MatLUFactor_SeqAIJ, 2927 0, 2928 MatSOR_SeqAIJ, 2929 MatTranspose_SeqAIJ, 2930 /*15*/ MatGetInfo_SeqAIJ, 2931 MatEqual_SeqAIJ, 2932 MatGetDiagonal_SeqAIJ, 2933 MatDiagonalScale_SeqAIJ, 2934 MatNorm_SeqAIJ, 2935 /*20*/ 0, 2936 MatAssemblyEnd_SeqAIJ, 2937 MatSetOption_SeqAIJ, 2938 MatZeroEntries_SeqAIJ, 2939 /*24*/ MatZeroRows_SeqAIJ, 2940 0, 2941 0, 2942 0, 2943 0, 2944 /*29*/ MatSetUp_SeqAIJ, 2945 0, 2946 0, 2947 MatGetArray_SeqAIJ, 2948 MatRestoreArray_SeqAIJ, 2949 /*34*/ MatDuplicate_SeqAIJ, 2950 0, 2951 0, 2952 MatILUFactor_SeqAIJ, 2953 0, 2954 /*39*/ MatAXPY_SeqAIJ, 2955 MatGetSubMatrices_SeqAIJ, 2956 MatIncreaseOverlap_SeqAIJ, 2957 MatGetValues_SeqAIJ, 2958 MatCopy_SeqAIJ, 2959 /*44*/ MatGetRowMax_SeqAIJ, 2960 MatScale_SeqAIJ, 2961 0, 2962 MatDiagonalSet_SeqAIJ, 2963 MatZeroRowsColumns_SeqAIJ, 2964 /*49*/ 0, 2965 MatGetRowIJ_SeqAIJ, 2966 MatRestoreRowIJ_SeqAIJ, 2967 MatGetColumnIJ_SeqAIJ, 2968 MatRestoreColumnIJ_SeqAIJ, 2969 /*54*/ MatFDColoringCreate_SeqAIJ, 2970 0, 2971 0, 2972 MatPermute_SeqAIJ, 2973 0, 2974 /*59*/ 0, 2975 MatDestroy_SeqAIJ, 2976 MatView_SeqAIJ, 2977 0, 2978 0, 2979 /*64*/ 0, 2980 0, 2981 0, 2982 0, 2983 0, 2984 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2985 MatGetRowMinAbs_SeqAIJ, 2986 0, 2987 MatSetColoring_SeqAIJ, 2988 #if defined(PETSC_HAVE_ADIC) 2989 MatSetValuesAdic_SeqAIJ, 2990 #else 2991 0, 2992 #endif 2993 /*74*/ MatSetValuesAdifor_SeqAIJ, 2994 MatFDColoringApply_AIJ, 2995 0, 2996 0, 2997 0, 2998 /*79*/ MatFindZeroDiagonals_SeqAIJ, 2999 0, 3000 0, 3001 0, 3002 MatLoad_SeqAIJ, 3003 /*84*/ MatIsSymmetric_SeqAIJ, 3004 MatIsHermitian_SeqAIJ, 3005 0, 3006 0, 3007 0, 3008 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 3009 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3010 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3011 MatPtAP_Basic, 3012 MatPtAPSymbolic_SeqAIJ, 3013 /*94*/ MatPtAPNumeric_SeqAIJ, 3014 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3015 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3016 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3017 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 3018 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3019 0, 3020 0, 3021 MatConjugate_SeqAIJ, 3022 0, 3023 /*104*/MatSetValuesRow_SeqAIJ, 3024 MatRealPart_SeqAIJ, 3025 MatImaginaryPart_SeqAIJ, 3026 0, 3027 0, 3028 /*109*/MatMatSolve_SeqAIJ, 3029 0, 3030 MatGetRowMin_SeqAIJ, 3031 0, 3032 MatMissingDiagonal_SeqAIJ, 3033 /*114*/0, 3034 0, 3035 0, 3036 0, 3037 0, 3038 /*119*/0, 3039 0, 3040 0, 3041 0, 3042 MatGetMultiProcBlock_SeqAIJ, 3043 /*124*/MatFindNonzeroRows_SeqAIJ, 3044 MatGetColumnNorms_SeqAIJ, 3045 MatInvertBlockDiagonal_SeqAIJ, 3046 0, 3047 0, 3048 /*129*/0, 3049 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3050 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3051 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3052 MatTransposeColoringCreate_SeqAIJ, 3053 /*134*/MatTransColoringApplySpToDen_SeqAIJ, 3054 MatTransColoringApplyDenToSp_SeqAIJ, 3055 MatRARt_SeqAIJ_SeqAIJ, 3056 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3057 MatRARtNumeric_SeqAIJ_SeqAIJ 3058 }; 3059 3060 EXTERN_C_BEGIN 3061 #undef __FUNCT__ 3062 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3063 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3064 { 3065 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3066 PetscInt i,nz,n; 3067 3068 PetscFunctionBegin; 3069 3070 nz = aij->maxnz; 3071 n = mat->rmap->n; 3072 for (i=0; i<nz; i++) { 3073 aij->j[i] = indices[i]; 3074 } 3075 aij->nz = nz; 3076 for (i=0; i<n; i++) { 3077 aij->ilen[i] = aij->imax[i]; 3078 } 3079 3080 PetscFunctionReturn(0); 3081 } 3082 EXTERN_C_END 3083 3084 #undef __FUNCT__ 3085 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3086 /*@ 3087 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3088 in the matrix. 3089 3090 Input Parameters: 3091 + mat - the SeqAIJ matrix 3092 - indices - the column indices 3093 3094 Level: advanced 3095 3096 Notes: 3097 This can be called if you have precomputed the nonzero structure of the 3098 matrix and want to provide it to the matrix object to improve the performance 3099 of the MatSetValues() operation. 3100 3101 You MUST have set the correct numbers of nonzeros per row in the call to 3102 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3103 3104 MUST be called before any calls to MatSetValues(); 3105 3106 The indices should start with zero, not one. 3107 3108 @*/ 3109 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3110 { 3111 PetscErrorCode ierr; 3112 3113 PetscFunctionBegin; 3114 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3115 PetscValidPointer(indices,2); 3116 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 3117 PetscFunctionReturn(0); 3118 } 3119 3120 /* ----------------------------------------------------------------------------------------*/ 3121 3122 EXTERN_C_BEGIN 3123 #undef __FUNCT__ 3124 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3125 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3126 { 3127 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3128 PetscErrorCode ierr; 3129 size_t nz = aij->i[mat->rmap->n]; 3130 3131 PetscFunctionBegin; 3132 if (aij->nonew != 1) { 3133 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3134 } 3135 3136 /* allocate space for values if not already there */ 3137 if (!aij->saved_values) { 3138 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 3139 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3140 } 3141 3142 /* copy values over */ 3143 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3144 PetscFunctionReturn(0); 3145 } 3146 EXTERN_C_END 3147 3148 #undef __FUNCT__ 3149 #define __FUNCT__ "MatStoreValues" 3150 /*@ 3151 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3152 example, reuse of the linear part of a Jacobian, while recomputing the 3153 nonlinear portion. 3154 3155 Collect on Mat 3156 3157 Input Parameters: 3158 . mat - the matrix (currently only AIJ matrices support this option) 3159 3160 Level: advanced 3161 3162 Common Usage, with SNESSolve(): 3163 $ Create Jacobian matrix 3164 $ Set linear terms into matrix 3165 $ Apply boundary conditions to matrix, at this time matrix must have 3166 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3167 $ boundary conditions again will not change the nonzero structure 3168 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3169 $ ierr = MatStoreValues(mat); 3170 $ Call SNESSetJacobian() with matrix 3171 $ In your Jacobian routine 3172 $ ierr = MatRetrieveValues(mat); 3173 $ Set nonlinear terms in matrix 3174 3175 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3176 $ // build linear portion of Jacobian 3177 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3178 $ ierr = MatStoreValues(mat); 3179 $ loop over nonlinear iterations 3180 $ ierr = MatRetrieveValues(mat); 3181 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3182 $ // call MatAssemblyBegin/End() on matrix 3183 $ Solve linear system with Jacobian 3184 $ endloop 3185 3186 Notes: 3187 Matrix must already be assemblied before calling this routine 3188 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3189 calling this routine. 3190 3191 When this is called multiple times it overwrites the previous set of stored values 3192 and does not allocated additional space. 3193 3194 .seealso: MatRetrieveValues() 3195 3196 @*/ 3197 PetscErrorCode MatStoreValues(Mat mat) 3198 { 3199 PetscErrorCode ierr; 3200 3201 PetscFunctionBegin; 3202 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3203 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3204 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3205 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3206 PetscFunctionReturn(0); 3207 } 3208 3209 EXTERN_C_BEGIN 3210 #undef __FUNCT__ 3211 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3212 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3213 { 3214 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3215 PetscErrorCode ierr; 3216 PetscInt nz = aij->i[mat->rmap->n]; 3217 3218 PetscFunctionBegin; 3219 if (aij->nonew != 1) { 3220 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3221 } 3222 if (!aij->saved_values) { 3223 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3224 } 3225 /* copy values over */ 3226 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3227 PetscFunctionReturn(0); 3228 } 3229 EXTERN_C_END 3230 3231 #undef __FUNCT__ 3232 #define __FUNCT__ "MatRetrieveValues" 3233 /*@ 3234 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3235 example, reuse of the linear part of a Jacobian, while recomputing the 3236 nonlinear portion. 3237 3238 Collect on Mat 3239 3240 Input Parameters: 3241 . mat - the matrix (currently on AIJ matrices support this option) 3242 3243 Level: advanced 3244 3245 .seealso: MatStoreValues() 3246 3247 @*/ 3248 PetscErrorCode MatRetrieveValues(Mat mat) 3249 { 3250 PetscErrorCode ierr; 3251 3252 PetscFunctionBegin; 3253 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3254 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3255 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3256 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3257 PetscFunctionReturn(0); 3258 } 3259 3260 3261 /* --------------------------------------------------------------------------------*/ 3262 #undef __FUNCT__ 3263 #define __FUNCT__ "MatCreateSeqAIJ" 3264 /*@C 3265 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3266 (the default parallel PETSc format). For good matrix assembly performance 3267 the user should preallocate the matrix storage by setting the parameter nz 3268 (or the array nnz). By setting these parameters accurately, performance 3269 during matrix assembly can be increased by more than a factor of 50. 3270 3271 Collective on MPI_Comm 3272 3273 Input Parameters: 3274 + comm - MPI communicator, set to PETSC_COMM_SELF 3275 . m - number of rows 3276 . n - number of columns 3277 . nz - number of nonzeros per row (same for all rows) 3278 - nnz - array containing the number of nonzeros in the various rows 3279 (possibly different for each row) or PETSC_NULL 3280 3281 Output Parameter: 3282 . A - the matrix 3283 3284 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3285 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3286 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3287 3288 Notes: 3289 If nnz is given then nz is ignored 3290 3291 The AIJ format (also called the Yale sparse matrix format or 3292 compressed row storage), is fully compatible with standard Fortran 77 3293 storage. That is, the stored row and column indices can begin at 3294 either one (as in Fortran) or zero. See the users' manual for details. 3295 3296 Specify the preallocated storage with either nz or nnz (not both). 3297 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3298 allocation. For large problems you MUST preallocate memory or you 3299 will get TERRIBLE performance, see the users' manual chapter on matrices. 3300 3301 By default, this format uses inodes (identical nodes) when possible, to 3302 improve numerical efficiency of matrix-vector products and solves. We 3303 search for consecutive rows with the same nonzero structure, thereby 3304 reusing matrix information to achieve increased efficiency. 3305 3306 Options Database Keys: 3307 + -mat_no_inode - Do not use inodes 3308 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3309 3310 Level: intermediate 3311 3312 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3313 3314 @*/ 3315 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3316 { 3317 PetscErrorCode ierr; 3318 3319 PetscFunctionBegin; 3320 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3321 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3322 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3323 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3324 PetscFunctionReturn(0); 3325 } 3326 3327 #undef __FUNCT__ 3328 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3329 /*@C 3330 MatSeqAIJSetPreallocation - For good matrix assembly performance 3331 the user should preallocate the matrix storage by setting the parameter nz 3332 (or the array nnz). By setting these parameters accurately, performance 3333 during matrix assembly can be increased by more than a factor of 50. 3334 3335 Collective on MPI_Comm 3336 3337 Input Parameters: 3338 + B - The matrix-free 3339 . nz - number of nonzeros per row (same for all rows) 3340 - nnz - array containing the number of nonzeros in the various rows 3341 (possibly different for each row) or PETSC_NULL 3342 3343 Notes: 3344 If nnz is given then nz is ignored 3345 3346 The AIJ format (also called the Yale sparse matrix format or 3347 compressed row storage), is fully compatible with standard Fortran 77 3348 storage. That is, the stored row and column indices can begin at 3349 either one (as in Fortran) or zero. See the users' manual for details. 3350 3351 Specify the preallocated storage with either nz or nnz (not both). 3352 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3353 allocation. For large problems you MUST preallocate memory or you 3354 will get TERRIBLE performance, see the users' manual chapter on matrices. 3355 3356 You can call MatGetInfo() to get information on how effective the preallocation was; 3357 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3358 You can also run with the option -info and look for messages with the string 3359 malloc in them to see if additional memory allocation was needed. 3360 3361 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3362 entries or columns indices 3363 3364 By default, this format uses inodes (identical nodes) when possible, to 3365 improve numerical efficiency of matrix-vector products and solves. We 3366 search for consecutive rows with the same nonzero structure, thereby 3367 reusing matrix information to achieve increased efficiency. 3368 3369 Options Database Keys: 3370 + -mat_no_inode - Do not use inodes 3371 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3372 - -mat_aij_oneindex - Internally use indexing starting at 1 3373 rather than 0. Note that when calling MatSetValues(), 3374 the user still MUST index entries starting at 0! 3375 3376 Level: intermediate 3377 3378 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3379 3380 @*/ 3381 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3382 { 3383 PetscErrorCode ierr; 3384 3385 PetscFunctionBegin; 3386 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3387 PetscValidType(B,1); 3388 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3389 PetscFunctionReturn(0); 3390 } 3391 3392 EXTERN_C_BEGIN 3393 #undef __FUNCT__ 3394 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3395 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3396 { 3397 Mat_SeqAIJ *b; 3398 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3399 PetscErrorCode ierr; 3400 PetscInt i; 3401 3402 PetscFunctionBegin; 3403 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3404 if (nz == MAT_SKIP_ALLOCATION) { 3405 skipallocation = PETSC_TRUE; 3406 nz = 0; 3407 } 3408 3409 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3410 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3411 3412 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3413 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3414 if (nnz) { 3415 for (i=0; i<B->rmap->n; i++) { 3416 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]); 3417 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); 3418 } 3419 } 3420 3421 B->preallocated = PETSC_TRUE; 3422 b = (Mat_SeqAIJ*)B->data; 3423 3424 if (!skipallocation) { 3425 if (!b->imax) { 3426 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3427 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3428 } 3429 if (!nnz) { 3430 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3431 else if (nz < 0) nz = 1; 3432 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3433 nz = nz*B->rmap->n; 3434 } else { 3435 nz = 0; 3436 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3437 } 3438 /* b->ilen will count nonzeros in each row so far. */ 3439 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3440 3441 /* allocate the matrix space */ 3442 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3443 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3444 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3445 b->i[0] = 0; 3446 for (i=1; i<B->rmap->n+1; i++) { 3447 b->i[i] = b->i[i-1] + b->imax[i-1]; 3448 } 3449 b->singlemalloc = PETSC_TRUE; 3450 b->free_a = PETSC_TRUE; 3451 b->free_ij = PETSC_TRUE; 3452 } else { 3453 b->free_a = PETSC_FALSE; 3454 b->free_ij = PETSC_FALSE; 3455 } 3456 3457 b->nz = 0; 3458 b->maxnz = nz; 3459 B->info.nz_unneeded = (double)b->maxnz; 3460 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 3461 PetscFunctionReturn(0); 3462 } 3463 EXTERN_C_END 3464 3465 #undef __FUNCT__ 3466 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3467 /*@ 3468 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3469 3470 Input Parameters: 3471 + B - the matrix 3472 . i - the indices into j for the start of each row (starts with zero) 3473 . j - the column indices for each row (starts with zero) these must be sorted for each row 3474 - v - optional values in the matrix 3475 3476 Level: developer 3477 3478 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3479 3480 .keywords: matrix, aij, compressed row, sparse, sequential 3481 3482 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3483 @*/ 3484 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3485 { 3486 PetscErrorCode ierr; 3487 3488 PetscFunctionBegin; 3489 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3490 PetscValidType(B,1); 3491 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3492 PetscFunctionReturn(0); 3493 } 3494 3495 EXTERN_C_BEGIN 3496 #undef __FUNCT__ 3497 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3498 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3499 { 3500 PetscInt i; 3501 PetscInt m,n; 3502 PetscInt nz; 3503 PetscInt *nnz, nz_max = 0; 3504 PetscScalar *values; 3505 PetscErrorCode ierr; 3506 3507 PetscFunctionBegin; 3508 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3509 3510 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3511 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3512 3513 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3514 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3515 for(i = 0; i < m; i++) { 3516 nz = Ii[i+1]- Ii[i]; 3517 nz_max = PetscMax(nz_max, nz); 3518 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3519 nnz[i] = nz; 3520 } 3521 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3522 ierr = PetscFree(nnz);CHKERRQ(ierr); 3523 3524 if (v) { 3525 values = (PetscScalar*) v; 3526 } else { 3527 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3528 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3529 } 3530 3531 for(i = 0; i < m; i++) { 3532 nz = Ii[i+1] - Ii[i]; 3533 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3534 } 3535 3536 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3537 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3538 3539 if (!v) { 3540 ierr = PetscFree(values);CHKERRQ(ierr); 3541 } 3542 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3543 PetscFunctionReturn(0); 3544 } 3545 EXTERN_C_END 3546 3547 #include <../src/mat/impls/dense/seq/dense.h> 3548 #include <petsc-private/petscaxpy.h> 3549 3550 #undef __FUNCT__ 3551 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3552 /* 3553 Computes (B'*A')' since computing B*A directly is untenable 3554 3555 n p p 3556 ( ) ( ) ( ) 3557 m ( A ) * n ( B ) = m ( C ) 3558 ( ) ( ) ( ) 3559 3560 */ 3561 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3562 { 3563 PetscErrorCode ierr; 3564 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3565 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3566 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3567 PetscInt i,n,m,q,p; 3568 const PetscInt *ii,*idx; 3569 const PetscScalar *b,*a,*a_q; 3570 PetscScalar *c,*c_q; 3571 3572 PetscFunctionBegin; 3573 m = A->rmap->n; 3574 n = A->cmap->n; 3575 p = B->cmap->n; 3576 a = sub_a->v; 3577 b = sub_b->a; 3578 c = sub_c->v; 3579 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3580 3581 ii = sub_b->i; 3582 idx = sub_b->j; 3583 for (i=0; i<n; i++) { 3584 q = ii[i+1] - ii[i]; 3585 while (q-->0) { 3586 c_q = c + m*(*idx); 3587 a_q = a + m*i; 3588 PetscAXPY(c_q,*b,a_q,m); 3589 idx++; 3590 b++; 3591 } 3592 } 3593 PetscFunctionReturn(0); 3594 } 3595 3596 #undef __FUNCT__ 3597 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3598 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3599 { 3600 PetscErrorCode ierr; 3601 PetscInt m=A->rmap->n,n=B->cmap->n; 3602 Mat Cmat; 3603 3604 PetscFunctionBegin; 3605 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); 3606 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3607 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3608 ierr = MatSetBlockSizes(Cmat,A->rmap->bs,B->cmap->bs);CHKERRQ(ierr); 3609 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3610 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3611 Cmat->assembled = PETSC_TRUE; 3612 Cmat->ops->matmult = MatMatMult_SeqDense_SeqAIJ; 3613 *C = Cmat; 3614 PetscFunctionReturn(0); 3615 } 3616 3617 /* ----------------------------------------------------------------*/ 3618 #undef __FUNCT__ 3619 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3620 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3621 { 3622 PetscErrorCode ierr; 3623 3624 PetscFunctionBegin; 3625 if (scall == MAT_INITIAL_MATRIX){ 3626 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3627 } 3628 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3629 PetscFunctionReturn(0); 3630 } 3631 3632 3633 /*MC 3634 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3635 based on compressed sparse row format. 3636 3637 Options Database Keys: 3638 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3639 3640 Level: beginner 3641 3642 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3643 M*/ 3644 3645 EXTERN_C_BEGIN 3646 #if defined(PETSC_HAVE_PASTIX) 3647 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3648 #endif 3649 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3650 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3651 #endif 3652 extern PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3653 extern PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3654 extern PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3655 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool *); 3656 #if defined(PETSC_HAVE_MUMPS) 3657 extern PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3658 #endif 3659 #if defined(PETSC_HAVE_SUPERLU) 3660 extern PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3661 #endif 3662 #if defined(PETSC_HAVE_SUPERLU_DIST) 3663 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3664 #endif 3665 #if defined(PETSC_HAVE_SPOOLES) 3666 extern PetscErrorCode MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3667 #endif 3668 #if defined(PETSC_HAVE_UMFPACK) 3669 extern PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3670 #endif 3671 #if defined(PETSC_HAVE_CHOLMOD) 3672 extern PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3673 #endif 3674 #if defined(PETSC_HAVE_LUSOL) 3675 extern PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3676 #endif 3677 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3678 extern PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3679 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3680 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3681 #endif 3682 EXTERN_C_END 3683 3684 3685 EXTERN_C_BEGIN 3686 #undef __FUNCT__ 3687 #define __FUNCT__ "MatCreate_SeqAIJ" 3688 PetscErrorCode MatCreate_SeqAIJ(Mat B) 3689 { 3690 Mat_SeqAIJ *b; 3691 PetscErrorCode ierr; 3692 PetscMPIInt size; 3693 3694 PetscFunctionBegin; 3695 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3696 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3697 3698 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3699 B->data = (void*)b; 3700 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3701 b->row = 0; 3702 b->col = 0; 3703 b->icol = 0; 3704 b->reallocs = 0; 3705 b->ignorezeroentries = PETSC_FALSE; 3706 b->roworiented = PETSC_TRUE; 3707 b->nonew = 0; 3708 b->diag = 0; 3709 b->solve_work = 0; 3710 B->spptr = 0; 3711 b->saved_values = 0; 3712 b->idiag = 0; 3713 b->mdiag = 0; 3714 b->ssor_work = 0; 3715 b->omega = 1.0; 3716 b->fshift = 0.0; 3717 b->idiagvalid = PETSC_FALSE; 3718 b->ibdiagvalid = PETSC_FALSE; 3719 b->keepnonzeropattern = PETSC_FALSE; 3720 b->xtoy = 0; 3721 b->XtoY = 0; 3722 B->same_nonzero = PETSC_FALSE; 3723 3724 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3725 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3726 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C","MatGetFactor_seqaij_matlab",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3727 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3728 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3729 #endif 3730 #if defined(PETSC_HAVE_PASTIX) 3731 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C","MatGetFactor_seqaij_pastix",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3732 #endif 3733 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C","MatGetFactor_seqaij_essl",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3735 #endif 3736 #if defined(PETSC_HAVE_SUPERLU) 3737 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C","MatGetFactor_seqaij_superlu",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3738 #endif 3739 #if defined(PETSC_HAVE_SUPERLU_DIST) 3740 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C","MatGetFactor_seqaij_superlu_dist",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3741 #endif 3742 #if defined(PETSC_HAVE_SPOOLES) 3743 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C","MatGetFactor_seqaij_spooles",MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3744 #endif 3745 #if defined(PETSC_HAVE_MUMPS) 3746 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C","MatGetFactor_aij_mumps",MatGetFactor_aij_mumps);CHKERRQ(ierr); 3747 #endif 3748 #if defined(PETSC_HAVE_UMFPACK) 3749 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C","MatGetFactor_seqaij_umfpack",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3750 #endif 3751 #if defined(PETSC_HAVE_CHOLMOD) 3752 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C","MatGetFactor_seqaij_cholmod",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3753 #endif 3754 #if defined(PETSC_HAVE_LUSOL) 3755 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C","MatGetFactor_seqaij_lusol",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3756 #endif 3757 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C","MatGetFactor_seqaij_petsc",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3758 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C","MatGetFactorAvailable_seqaij_petsc",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3759 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C","MatGetFactor_seqaij_bas",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 3760 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C","MatSeqAIJSetColumnIndices_SeqAIJ",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3761 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C","MatStoreValues_SeqAIJ",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3762 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C","MatRetrieveValues_SeqAIJ",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3763 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C","MatConvert_SeqAIJ_SeqSBAIJ",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3764 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C","MatConvert_SeqAIJ_SeqBAIJ",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3765 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijperm_C","MatConvert_SeqAIJ_SeqAIJPERM",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3766 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C","MatConvert_SeqAIJ_SeqAIJCRL",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3767 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C","MatIsTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3768 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C","MatIsHermitianTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3769 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C","MatSeqAIJSetPreallocation_SeqAIJ",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3770 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C","MatSeqAIJSetPreallocationCSR_SeqAIJ",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3771 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C","MatReorderForNonzeroDiagonal_SeqAIJ",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3772 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C","MatMatMult_SeqDense_SeqAIJ",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3773 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C","MatMatMultSymbolic_SeqDense_SeqAIJ",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3774 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C","MatMatMultNumeric_SeqDense_SeqAIJ",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3775 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3776 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3777 PetscFunctionReturn(0); 3778 } 3779 EXTERN_C_END 3780 3781 #undef __FUNCT__ 3782 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3783 /* 3784 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3785 */ 3786 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3787 { 3788 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3789 PetscErrorCode ierr; 3790 PetscInt i,m = A->rmap->n; 3791 3792 PetscFunctionBegin; 3793 c = (Mat_SeqAIJ*)C->data; 3794 3795 C->factortype = A->factortype; 3796 c->row = 0; 3797 c->col = 0; 3798 c->icol = 0; 3799 c->reallocs = 0; 3800 3801 C->assembled = PETSC_TRUE; 3802 3803 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 3804 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 3805 3806 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3807 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3808 for (i=0; i<m; i++) { 3809 c->imax[i] = a->imax[i]; 3810 c->ilen[i] = a->ilen[i]; 3811 } 3812 3813 /* allocate the matrix space */ 3814 if (mallocmatspace){ 3815 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3816 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3817 c->singlemalloc = PETSC_TRUE; 3818 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3819 if (m > 0) { 3820 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3821 if (cpvalues == MAT_COPY_VALUES) { 3822 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3823 } else { 3824 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3825 } 3826 } 3827 } 3828 3829 c->ignorezeroentries = a->ignorezeroentries; 3830 c->roworiented = a->roworiented; 3831 c->nonew = a->nonew; 3832 if (a->diag) { 3833 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3834 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3835 for (i=0; i<m; i++) { 3836 c->diag[i] = a->diag[i]; 3837 } 3838 } else c->diag = 0; 3839 c->solve_work = 0; 3840 c->saved_values = 0; 3841 c->idiag = 0; 3842 c->ssor_work = 0; 3843 c->keepnonzeropattern = a->keepnonzeropattern; 3844 c->free_a = PETSC_TRUE; 3845 c->free_ij = PETSC_TRUE; 3846 c->xtoy = 0; 3847 c->XtoY = 0; 3848 3849 c->rmax = a->rmax; 3850 c->nz = a->nz; 3851 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3852 C->preallocated = PETSC_TRUE; 3853 3854 c->compressedrow.use = a->compressedrow.use; 3855 c->compressedrow.nrows = a->compressedrow.nrows; 3856 c->compressedrow.check = a->compressedrow.check; 3857 if (a->compressedrow.use){ 3858 i = a->compressedrow.nrows; 3859 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3860 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3861 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3862 } else { 3863 c->compressedrow.use = PETSC_FALSE; 3864 c->compressedrow.i = PETSC_NULL; 3865 c->compressedrow.rindex = PETSC_NULL; 3866 } 3867 C->same_nonzero = A->same_nonzero; 3868 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3869 3870 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3871 PetscFunctionReturn(0); 3872 } 3873 3874 #undef __FUNCT__ 3875 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3876 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3877 { 3878 PetscErrorCode ierr; 3879 3880 PetscFunctionBegin; 3881 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3882 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3883 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 3884 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 3885 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3886 PetscFunctionReturn(0); 3887 } 3888 3889 #undef __FUNCT__ 3890 #define __FUNCT__ "MatLoad_SeqAIJ" 3891 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 3892 { 3893 Mat_SeqAIJ *a; 3894 PetscErrorCode ierr; 3895 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3896 int fd; 3897 PetscMPIInt size; 3898 MPI_Comm comm; 3899 PetscInt bs = 1; 3900 3901 PetscFunctionBegin; 3902 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3903 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3904 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3905 3906 ierr = PetscOptionsBegin(comm,PETSC_NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 3907 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,PETSC_NULL);CHKERRQ(ierr); 3908 ierr = PetscOptionsEnd();CHKERRQ(ierr); 3909 3910 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3911 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3912 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3913 M = header[1]; N = header[2]; nz = header[3]; 3914 3915 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3916 3917 /* read in row lengths */ 3918 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3919 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3920 3921 /* check if sum of rowlengths is same as nz */ 3922 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3923 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); 3924 3925 /* set global size if not set already*/ 3926 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3927 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3928 } else { 3929 /* if sizes and type are already set, check if the vector global sizes are correct */ 3930 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3931 if (rows < 0 && cols < 0){ /* user might provide local size instead of global size */ 3932 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 3933 } 3934 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); 3935 } 3936 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3937 a = (Mat_SeqAIJ*)newMat->data; 3938 3939 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3940 3941 /* read in nonzero values */ 3942 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3943 3944 /* set matrix "i" values */ 3945 a->i[0] = 0; 3946 for (i=1; i<= M; i++) { 3947 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3948 a->ilen[i-1] = rowlengths[i-1]; 3949 } 3950 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3951 3952 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3953 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3954 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 3955 PetscFunctionReturn(0); 3956 } 3957 3958 #undef __FUNCT__ 3959 #define __FUNCT__ "MatEqual_SeqAIJ" 3960 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 3961 { 3962 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3963 PetscErrorCode ierr; 3964 #if defined(PETSC_USE_COMPLEX) 3965 PetscInt k; 3966 #endif 3967 3968 PetscFunctionBegin; 3969 /* If the matrix dimensions are not equal,or no of nonzeros */ 3970 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3971 *flg = PETSC_FALSE; 3972 PetscFunctionReturn(0); 3973 } 3974 3975 /* if the a->i are the same */ 3976 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3977 if (!*flg) PetscFunctionReturn(0); 3978 3979 /* if a->j are the same */ 3980 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3981 if (!*flg) PetscFunctionReturn(0); 3982 3983 /* if a->a are the same */ 3984 #if defined(PETSC_USE_COMPLEX) 3985 for (k=0; k<a->nz; k++){ 3986 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3987 *flg = PETSC_FALSE; 3988 PetscFunctionReturn(0); 3989 } 3990 } 3991 #else 3992 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3993 #endif 3994 PetscFunctionReturn(0); 3995 } 3996 3997 #undef __FUNCT__ 3998 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3999 /*@ 4000 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4001 provided by the user. 4002 4003 Collective on MPI_Comm 4004 4005 Input Parameters: 4006 + comm - must be an MPI communicator of size 1 4007 . m - number of rows 4008 . n - number of columns 4009 . i - row indices 4010 . j - column indices 4011 - a - matrix values 4012 4013 Output Parameter: 4014 . mat - the matrix 4015 4016 Level: intermediate 4017 4018 Notes: 4019 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4020 once the matrix is destroyed and not before 4021 4022 You cannot set new nonzero locations into this matrix, that will generate an error. 4023 4024 The i and j indices are 0 based 4025 4026 The format which is used for the sparse matrix input, is equivalent to a 4027 row-major ordering.. i.e for the following matrix, the input data expected is 4028 as shown: 4029 4030 1 0 0 4031 2 0 3 4032 4 5 6 4033 4034 i = {0,1,3,6} [size = nrow+1 = 3+1] 4035 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4036 v = {1,2,3,4,5,6} [size = nz = 6] 4037 4038 4039 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4040 4041 @*/ 4042 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 4043 { 4044 PetscErrorCode ierr; 4045 PetscInt ii; 4046 Mat_SeqAIJ *aij; 4047 #if defined(PETSC_USE_DEBUG) 4048 PetscInt jj; 4049 #endif 4050 4051 PetscFunctionBegin; 4052 if (i[0]) { 4053 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4054 } 4055 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4056 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4057 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4058 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4059 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4060 aij = (Mat_SeqAIJ*)(*mat)->data; 4061 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 4062 4063 aij->i = i; 4064 aij->j = j; 4065 aij->a = a; 4066 aij->singlemalloc = PETSC_FALSE; 4067 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4068 aij->free_a = PETSC_FALSE; 4069 aij->free_ij = PETSC_FALSE; 4070 4071 for (ii=0; ii<m; ii++) { 4072 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4073 #if defined(PETSC_USE_DEBUG) 4074 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]); 4075 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4076 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); 4077 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); 4078 } 4079 #endif 4080 } 4081 #if defined(PETSC_USE_DEBUG) 4082 for (ii=0; ii<aij->i[m]; ii++) { 4083 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4084 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]); 4085 } 4086 #endif 4087 4088 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4089 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4090 PetscFunctionReturn(0); 4091 } 4092 #undef __FUNCT__ 4093 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4094 /*@C 4095 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4096 provided by the user. 4097 4098 Collective on MPI_Comm 4099 4100 Input Parameters: 4101 + comm - must be an MPI communicator of size 1 4102 . m - number of rows 4103 . n - number of columns 4104 . i - row indices 4105 . j - column indices 4106 . a - matrix values 4107 . nz - number of nonzeros 4108 - idx - 0 or 1 based 4109 4110 Output Parameter: 4111 . mat - the matrix 4112 4113 Level: intermediate 4114 4115 Notes: 4116 The i and j indices are 0 based 4117 4118 The format which is used for the sparse matrix input, is equivalent to a 4119 row-major ordering.. i.e for the following matrix, the input data expected is 4120 as shown: 4121 4122 1 0 0 4123 2 0 3 4124 4 5 6 4125 4126 i = {0,1,1,2,2,2} 4127 j = {0,0,2,0,1,2} 4128 v = {1,2,3,4,5,6} 4129 4130 4131 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4132 4133 @*/ 4134 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4135 { 4136 PetscErrorCode ierr; 4137 PetscInt ii, *nnz, one = 1,row,col; 4138 4139 4140 PetscFunctionBegin; 4141 ierr = PetscMalloc(m*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 4142 ierr = PetscMemzero(nnz,m*sizeof(PetscInt));CHKERRQ(ierr); 4143 for (ii = 0; ii < nz; ii++){ 4144 nnz[i[ii]] += 1; 4145 } 4146 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4147 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4148 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4149 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4150 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4151 for (ii = 0; ii < nz; ii++){ 4152 if (idx){ 4153 row = i[ii] - 1; 4154 col = j[ii] - 1; 4155 } else { 4156 row = i[ii]; 4157 col = j[ii]; 4158 } 4159 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4160 } 4161 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4162 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4163 ierr = PetscFree(nnz);CHKERRQ(ierr); 4164 PetscFunctionReturn(0); 4165 } 4166 4167 #undef __FUNCT__ 4168 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4169 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4170 { 4171 PetscErrorCode ierr; 4172 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4173 4174 PetscFunctionBegin; 4175 if (coloring->ctype == IS_COLORING_GLOBAL) { 4176 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4177 a->coloring = coloring; 4178 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4179 PetscInt i,*larray; 4180 ISColoring ocoloring; 4181 ISColoringValue *colors; 4182 4183 /* set coloring for diagonal portion */ 4184 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 4185 for (i=0; i<A->cmap->n; i++) { 4186 larray[i] = i; 4187 } 4188 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 4189 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 4190 for (i=0; i<A->cmap->n; i++) { 4191 colors[i] = coloring->colors[larray[i]]; 4192 } 4193 ierr = PetscFree(larray);CHKERRQ(ierr); 4194 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4195 a->coloring = ocoloring; 4196 } 4197 PetscFunctionReturn(0); 4198 } 4199 4200 #if defined(PETSC_HAVE_ADIC) 4201 EXTERN_C_BEGIN 4202 #include <adic/ad_utils.h> 4203 EXTERN_C_END 4204 4205 #undef __FUNCT__ 4206 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 4207 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 4208 { 4209 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4210 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 4211 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 4212 ISColoringValue *color; 4213 4214 PetscFunctionBegin; 4215 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4216 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 4217 color = a->coloring->colors; 4218 /* loop over rows */ 4219 for (i=0; i<m; i++) { 4220 nz = ii[i+1] - ii[i]; 4221 /* loop over columns putting computed value into matrix */ 4222 for (j=0; j<nz; j++) { 4223 *v++ = values[color[*jj++]]; 4224 } 4225 values += nlen; /* jump to next row of derivatives */ 4226 } 4227 PetscFunctionReturn(0); 4228 } 4229 #endif 4230 4231 #undef __FUNCT__ 4232 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4233 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4234 { 4235 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4236 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4237 MatScalar *v = a->a; 4238 PetscScalar *values = (PetscScalar *)advalues; 4239 ISColoringValue *color; 4240 4241 PetscFunctionBegin; 4242 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4243 color = a->coloring->colors; 4244 /* loop over rows */ 4245 for (i=0; i<m; i++) { 4246 nz = ii[i+1] - ii[i]; 4247 /* loop over columns putting computed value into matrix */ 4248 for (j=0; j<nz; j++) { 4249 *v++ = values[color[*jj++]]; 4250 } 4251 values += nl; /* jump to next row of derivatives */ 4252 } 4253 PetscFunctionReturn(0); 4254 } 4255 4256 /* 4257 Special version for direct calls from Fortran 4258 */ 4259 #include <petsc-private/fortranimpl.h> 4260 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4261 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4262 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4263 #define matsetvaluesseqaij_ matsetvaluesseqaij 4264 #endif 4265 4266 /* Change these macros so can be used in void function */ 4267 #undef CHKERRQ 4268 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 4269 #undef SETERRQ2 4270 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4271 #undef SETERRQ3 4272 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4273 4274 EXTERN_C_BEGIN 4275 #undef __FUNCT__ 4276 #define __FUNCT__ "matsetvaluesseqaij_" 4277 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 4278 { 4279 Mat A = *AA; 4280 PetscInt m = *mm, n = *nn; 4281 InsertMode is = *isis; 4282 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4283 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4284 PetscInt *imax,*ai,*ailen; 4285 PetscErrorCode ierr; 4286 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4287 MatScalar *ap,value,*aa; 4288 PetscBool ignorezeroentries = a->ignorezeroentries; 4289 PetscBool roworiented = a->roworiented; 4290 4291 PetscFunctionBegin; 4292 MatCheckPreallocated(A,1); 4293 imax = a->imax; 4294 ai = a->i; 4295 ailen = a->ilen; 4296 aj = a->j; 4297 aa = a->a; 4298 4299 for (k=0; k<m; k++) { /* loop over added rows */ 4300 row = im[k]; 4301 if (row < 0) continue; 4302 #if defined(PETSC_USE_DEBUG) 4303 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4304 #endif 4305 rp = aj + ai[row]; ap = aa + ai[row]; 4306 rmax = imax[row]; nrow = ailen[row]; 4307 low = 0; 4308 high = nrow; 4309 for (l=0; l<n; l++) { /* loop over added columns */ 4310 if (in[l] < 0) continue; 4311 #if defined(PETSC_USE_DEBUG) 4312 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4313 #endif 4314 col = in[l]; 4315 if (roworiented) { 4316 value = v[l + k*n]; 4317 } else { 4318 value = v[k + l*m]; 4319 } 4320 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4321 4322 if (col <= lastcol) low = 0; else high = nrow; 4323 lastcol = col; 4324 while (high-low > 5) { 4325 t = (low+high)/2; 4326 if (rp[t] > col) high = t; 4327 else low = t; 4328 } 4329 for (i=low; i<high; i++) { 4330 if (rp[i] > col) break; 4331 if (rp[i] == col) { 4332 if (is == ADD_VALUES) ap[i] += value; 4333 else ap[i] = value; 4334 goto noinsert; 4335 } 4336 } 4337 if (value == 0.0 && ignorezeroentries) goto noinsert; 4338 if (nonew == 1) goto noinsert; 4339 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4340 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4341 N = nrow++ - 1; a->nz++; high++; 4342 /* shift up all the later entries in this row */ 4343 for (ii=N; ii>=i; ii--) { 4344 rp[ii+1] = rp[ii]; 4345 ap[ii+1] = ap[ii]; 4346 } 4347 rp[i] = col; 4348 ap[i] = value; 4349 noinsert:; 4350 low = i + 1; 4351 } 4352 ailen[row] = nrow; 4353 } 4354 A->same_nonzero = PETSC_FALSE; 4355 PetscFunctionReturnVoid(); 4356 } 4357 EXTERN_C_END 4358 4359