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