1 /*$Id: aij.c,v 1.385 2001/09/07 20:09:22 bsmith Exp $*/ 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 8 #include "src/vec/vecimpl.h" 9 #include "src/inline/spops.h" 10 #include "src/inline/dot.h" 11 #include "petscbt.h" 12 13 14 EXTERN int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 18 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int *ia[],int *ja[],PetscTruth *done) 19 { 20 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 21 int ierr,i,ishift; 22 23 PetscFunctionBegin; 24 *m = A->m; 25 if (!ia) PetscFunctionReturn(0); 26 ishift = 0; 27 if (symmetric && !A->structurally_symmetric) { 28 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 29 } else if (oshift == 1) { 30 int nz = a->i[A->m]; 31 /* malloc space and add 1 to i and j indices */ 32 ierr = PetscMalloc((A->m+1)*sizeof(int),ia);CHKERRQ(ierr); 33 ierr = PetscMalloc((nz+1)*sizeof(int),ja);CHKERRQ(ierr); 34 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 35 for (i=0; i<A->m+1; i++) (*ia)[i] = a->i[i] + 1; 36 } else { 37 *ia = a->i; *ja = a->j; 38 } 39 PetscFunctionReturn(0); 40 } 41 42 #undef __FUNCT__ 43 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 44 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int *ia[],int *ja[],PetscTruth *done) 45 { 46 int ierr; 47 48 PetscFunctionBegin; 49 if (!ia) PetscFunctionReturn(0); 50 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 51 ierr = PetscFree(*ia);CHKERRQ(ierr); 52 ierr = PetscFree(*ja);CHKERRQ(ierr); 53 } 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 59 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int *ia[],int *ja[],PetscTruth *done) 60 { 61 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 62 int ierr,i,*collengths,*cia,*cja,n = A->n,m = A->m; 63 int nz = a->i[m],row,*jj,mr,col; 64 65 PetscFunctionBegin; 66 *nn = A->n; 67 if (!ia) PetscFunctionReturn(0); 68 if (symmetric) { 69 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 70 } else { 71 ierr = PetscMalloc((n+1)*sizeof(int),&collengths);CHKERRQ(ierr); 72 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 73 ierr = PetscMalloc((n+1)*sizeof(int),&cia);CHKERRQ(ierr); 74 ierr = PetscMalloc((nz+1)*sizeof(int),&cja);CHKERRQ(ierr); 75 jj = a->j; 76 for (i=0; i<nz; i++) { 77 collengths[jj[i]]++; 78 } 79 cia[0] = oshift; 80 for (i=0; i<n; i++) { 81 cia[i+1] = cia[i] + collengths[i]; 82 } 83 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 84 jj = a->j; 85 for (row=0; row<m; row++) { 86 mr = a->i[row+1] - a->i[row]; 87 for (i=0; i<mr; i++) { 88 col = *jj++; 89 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 90 } 91 } 92 ierr = PetscFree(collengths);CHKERRQ(ierr); 93 *ia = cia; *ja = cja; 94 } 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 100 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int *ia[],int *ja[],PetscTruth *done) 101 { 102 int ierr; 103 104 PetscFunctionBegin; 105 if (!ia) PetscFunctionReturn(0); 106 107 ierr = PetscFree(*ia);CHKERRQ(ierr); 108 ierr = PetscFree(*ja);CHKERRQ(ierr); 109 110 PetscFunctionReturn(0); 111 } 112 113 #define CHUNKSIZE 15 114 115 #undef __FUNCT__ 116 #define __FUNCT__ "MatSetValues_SeqAIJ" 117 int MatSetValues_SeqAIJ(Mat A,int m,const int im[],int n,const int in[],const PetscScalar v[],InsertMode is) 118 { 119 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 120 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted = a->sorted; 121 int *imax = a->imax,*ai = a->i,*ailen = a->ilen; 122 int *aj = a->j,nonew = a->nonew,ierr; 123 PetscScalar *ap,value,*aa = a->a; 124 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 125 PetscTruth roworiented = a->roworiented; 126 127 PetscFunctionBegin; 128 for (k=0; k<m; k++) { /* loop over added rows */ 129 row = im[k]; 130 if (row < 0) continue; 131 #if defined(PETSC_USE_BOPT_g) 132 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m-1); 133 #endif 134 rp = aj + ai[row]; ap = aa + ai[row]; 135 rmax = imax[row]; nrow = ailen[row]; 136 low = 0; 137 for (l=0; l<n; l++) { /* loop over added columns */ 138 if (in[l] < 0) continue; 139 #if defined(PETSC_USE_BOPT_g) 140 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->n-1); 141 #endif 142 col = in[l]; 143 if (roworiented) { 144 value = v[l + k*n]; 145 } else { 146 value = v[k + l*m]; 147 } 148 if (value == 0.0 && ignorezeroentries) continue; 149 150 if (!sorted) low = 0; high = nrow; 151 while (high-low > 5) { 152 t = (low+high)/2; 153 if (rp[t] > col) high = t; 154 else low = t; 155 } 156 for (i=low; i<high; i++) { 157 if (rp[i] > col) break; 158 if (rp[i] == col) { 159 if (is == ADD_VALUES) ap[i] += value; 160 else ap[i] = value; 161 goto noinsert; 162 } 163 } 164 if (nonew == 1) goto noinsert; 165 else if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%d,%d) in the matrix",row,col); 166 if (nrow >= rmax) { 167 /* there is no extra room in row, therefore enlarge */ 168 int new_nz = ai[A->m] + CHUNKSIZE,*new_i,*new_j; 169 size_t len; 170 PetscScalar *new_a; 171 172 if (nonew == -2) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%d,%d) in the matrix requiring new malloc()",row,col); 173 174 /* malloc new storage space */ 175 len = ((size_t) new_nz)*(sizeof(int)+sizeof(PetscScalar))+(A->m+1)*sizeof(int); 176 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 177 new_j = (int*)(new_a + new_nz); 178 new_i = new_j + new_nz; 179 180 /* copy over old data into new slots */ 181 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 182 for (ii=row+1; ii<A->m+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 183 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int));CHKERRQ(ierr); 184 len = (((size_t) new_nz) - CHUNKSIZE - ai[row] - nrow ); 185 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int));CHKERRQ(ierr); 186 ierr = PetscMemcpy(new_a,aa,(((size_t) ai[row])+nrow)*sizeof(PetscScalar));CHKERRQ(ierr); 187 ierr = PetscMemcpy(new_a+ai[row]+nrow+CHUNKSIZE,aa+ai[row]+nrow,len*sizeof(PetscScalar));CHKERRQ(ierr); 188 /* free up old matrix storage */ 189 ierr = PetscFree(a->a);CHKERRQ(ierr); 190 if (!a->singlemalloc) { 191 ierr = PetscFree(a->i);CHKERRQ(ierr); 192 ierr = PetscFree(a->j);CHKERRQ(ierr); 193 } 194 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 195 a->singlemalloc = PETSC_TRUE; 196 197 rp = aj + ai[row]; ap = aa + ai[row] ; 198 rmax = imax[row] = imax[row] + CHUNKSIZE; 199 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(PetscScalar))); 200 a->maxnz += CHUNKSIZE; 201 a->reallocs++; 202 } 203 N = nrow++ - 1; a->nz++; 204 /* shift up all the later entries in this row */ 205 for (ii=N; ii>=i; ii--) { 206 rp[ii+1] = rp[ii]; 207 ap[ii+1] = ap[ii]; 208 } 209 rp[i] = col; 210 ap[i] = value; 211 noinsert:; 212 low = i + 1; 213 } 214 ailen[row] = nrow; 215 } 216 PetscFunctionReturn(0); 217 } 218 219 #undef __FUNCT__ 220 #define __FUNCT__ "MatGetValues_SeqAIJ" 221 int MatGetValues_SeqAIJ(Mat A,int m,const int im[],int n,const int in[],PetscScalar v[]) 222 { 223 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 224 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 225 int *ai = a->i,*ailen = a->ilen; 226 PetscScalar *ap,*aa = a->a,zero = 0.0; 227 228 PetscFunctionBegin; 229 for (k=0; k<m; k++) { /* loop over rows */ 230 row = im[k]; 231 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %d",row); 232 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m-1); 233 rp = aj + ai[row]; ap = aa + ai[row]; 234 nrow = ailen[row]; 235 for (l=0; l<n; l++) { /* loop over columns */ 236 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %d",in[l]); 237 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->n-1); 238 col = in[l] ; 239 high = nrow; low = 0; /* assume unsorted */ 240 while (high-low > 5) { 241 t = (low+high)/2; 242 if (rp[t] > col) high = t; 243 else low = t; 244 } 245 for (i=low; i<high; i++) { 246 if (rp[i] > col) break; 247 if (rp[i] == col) { 248 *v++ = ap[i]; 249 goto finished; 250 } 251 } 252 *v++ = zero; 253 finished:; 254 } 255 } 256 PetscFunctionReturn(0); 257 } 258 259 260 #undef __FUNCT__ 261 #define __FUNCT__ "MatView_SeqAIJ_Binary" 262 int MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 263 { 264 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 265 int i,fd,*col_lens,ierr; 266 267 PetscFunctionBegin; 268 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 269 ierr = PetscMalloc((4+A->m)*sizeof(int),&col_lens);CHKERRQ(ierr); 270 col_lens[0] = MAT_FILE_COOKIE; 271 col_lens[1] = A->m; 272 col_lens[2] = A->n; 273 col_lens[3] = a->nz; 274 275 /* store lengths of each row and write (including header) to file */ 276 for (i=0; i<A->m; i++) { 277 col_lens[4+i] = a->i[i+1] - a->i[i]; 278 } 279 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,1);CHKERRQ(ierr); 280 ierr = PetscFree(col_lens);CHKERRQ(ierr); 281 282 /* store column indices (zero start index) */ 283 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,0);CHKERRQ(ierr); 284 285 /* store nonzero values */ 286 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,0);CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 extern int MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 291 292 #undef __FUNCT__ 293 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 294 int MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 295 { 296 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 297 int ierr,i,j,m = A->m,shift=0; 298 char *name; 299 PetscViewerFormat format; 300 301 PetscFunctionBegin; 302 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 303 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 304 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_INFO) { 305 if (a->inode.size) { 306 ierr = PetscViewerASCIIPrintf(viewer,"using I-node routines: found %d nodes, limit used is %d\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr); 307 } else { 308 ierr = PetscViewerASCIIPrintf(viewer,"not using I-node routines\n");CHKERRQ(ierr); 309 } 310 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 311 int nofinalvalue = 0; 312 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->n-!shift)) { 313 nofinalvalue = 1; 314 } 315 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 316 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %d %d \n",m,A->n);CHKERRQ(ierr); 317 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %d \n",a->nz);CHKERRQ(ierr); 318 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%d,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 319 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 320 321 for (i=0; i<m; i++) { 322 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 323 #if defined(PETSC_USE_COMPLEX) 324 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); 325 #else 326 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 327 #endif 328 } 329 } 330 if (nofinalvalue) { 331 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e\n",m,A->n,0.0);CHKERRQ(ierr); 332 } 333 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 334 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 335 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 336 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_SINGLE) && !defined(PETSC_USE_COMPLEX) 337 ierr = MatSeqAIJFactorInfo_Matlab(A,viewer);CHKERRQ(ierr); 338 #endif 339 PetscFunctionReturn(0); 340 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 341 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 342 for (i=0; i<m; i++) { 343 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 344 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 345 #if defined(PETSC_USE_COMPLEX) 346 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 347 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 348 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 349 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 350 } else if (PetscRealPart(a->a[j]) != 0.0) { 351 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 352 } 353 #else 354 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 355 #endif 356 } 357 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 358 } 359 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 360 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 361 int nzd=0,fshift=1,*sptr; 362 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 363 ierr = PetscMalloc((m+1)*sizeof(int),&sptr);CHKERRQ(ierr); 364 for (i=0; i<m; i++) { 365 sptr[i] = nzd+1; 366 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 367 if (a->j[j] >= i) { 368 #if defined(PETSC_USE_COMPLEX) 369 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 370 #else 371 if (a->a[j] != 0.0) nzd++; 372 #endif 373 } 374 } 375 } 376 sptr[m] = nzd+1; 377 ierr = PetscViewerASCIIPrintf(viewer," %d %d\n\n",m,nzd);CHKERRQ(ierr); 378 for (i=0; i<m+1; i+=6) { 379 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);} 380 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);} 381 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);} 382 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 383 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 384 else {ierr = PetscViewerASCIIPrintf(viewer," %d\n",sptr[i]);CHKERRQ(ierr);} 385 } 386 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 387 ierr = PetscFree(sptr);CHKERRQ(ierr); 388 for (i=0; i<m; i++) { 389 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 390 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %d ",a->j[j]+fshift);CHKERRQ(ierr);} 391 } 392 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 393 } 394 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 395 for (i=0; i<m; i++) { 396 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 397 if (a->j[j] >= i) { 398 #if defined(PETSC_USE_COMPLEX) 399 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 400 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 401 } 402 #else 403 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 404 #endif 405 } 406 } 407 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 408 } 409 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 410 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 411 int cnt = 0,jcnt; 412 PetscScalar value; 413 414 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 415 for (i=0; i<m; i++) { 416 jcnt = 0; 417 for (j=0; j<A->n; j++) { 418 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 419 value = a->a[cnt++]; 420 jcnt++; 421 } else { 422 value = 0.0; 423 } 424 #if defined(PETSC_USE_COMPLEX) 425 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 426 #else 427 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 428 #endif 429 } 430 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 431 } 432 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 433 } else { 434 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 435 for (i=0; i<m; i++) { 436 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 437 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 438 #if defined(PETSC_USE_COMPLEX) 439 if (PetscImaginaryPart(a->a[j]) > 0.0) { 440 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 441 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 442 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 443 } else { 444 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 445 } 446 #else 447 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 448 #endif 449 } 450 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 451 } 452 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 453 } 454 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 455 PetscFunctionReturn(0); 456 } 457 458 #undef __FUNCT__ 459 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 460 int MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 461 { 462 Mat A = (Mat) Aa; 463 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 464 int ierr,i,j,m = A->m,color; 465 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 466 PetscViewer viewer; 467 PetscViewerFormat format; 468 469 PetscFunctionBegin; 470 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 471 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 472 473 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 474 /* loop over matrix elements drawing boxes */ 475 476 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 477 /* Blue for negative, Cyan for zero and Red for positive */ 478 color = PETSC_DRAW_BLUE; 479 for (i=0; i<m; i++) { 480 y_l = m - i - 1.0; y_r = y_l + 1.0; 481 for (j=a->i[i]; j<a->i[i+1]; j++) { 482 x_l = a->j[j] ; x_r = x_l + 1.0; 483 #if defined(PETSC_USE_COMPLEX) 484 if (PetscRealPart(a->a[j]) >= 0.) continue; 485 #else 486 if (a->a[j] >= 0.) continue; 487 #endif 488 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 489 } 490 } 491 color = PETSC_DRAW_CYAN; 492 for (i=0; i<m; i++) { 493 y_l = m - i - 1.0; y_r = y_l + 1.0; 494 for (j=a->i[i]; j<a->i[i+1]; j++) { 495 x_l = a->j[j]; x_r = x_l + 1.0; 496 if (a->a[j] != 0.) continue; 497 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 498 } 499 } 500 color = PETSC_DRAW_RED; 501 for (i=0; i<m; i++) { 502 y_l = m - i - 1.0; y_r = y_l + 1.0; 503 for (j=a->i[i]; j<a->i[i+1]; j++) { 504 x_l = a->j[j]; x_r = x_l + 1.0; 505 #if defined(PETSC_USE_COMPLEX) 506 if (PetscRealPart(a->a[j]) <= 0.) continue; 507 #else 508 if (a->a[j] <= 0.) continue; 509 #endif 510 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 511 } 512 } 513 } else { 514 /* use contour shading to indicate magnitude of values */ 515 /* first determine max of all nonzero values */ 516 int nz = a->nz,count; 517 PetscDraw popup; 518 PetscReal scale; 519 520 for (i=0; i<nz; i++) { 521 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 522 } 523 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 524 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 525 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 526 count = 0; 527 for (i=0; i<m; i++) { 528 y_l = m - i - 1.0; y_r = y_l + 1.0; 529 for (j=a->i[i]; j<a->i[i+1]; j++) { 530 x_l = a->j[j]; x_r = x_l + 1.0; 531 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 532 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 533 count++; 534 } 535 } 536 } 537 PetscFunctionReturn(0); 538 } 539 540 #undef __FUNCT__ 541 #define __FUNCT__ "MatView_SeqAIJ_Draw" 542 int MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 543 { 544 int ierr; 545 PetscDraw draw; 546 PetscReal xr,yr,xl,yl,h,w; 547 PetscTruth isnull; 548 549 PetscFunctionBegin; 550 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 551 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 552 if (isnull) PetscFunctionReturn(0); 553 554 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 555 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 556 xr += w; yr += h; xl = -w; yl = -h; 557 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 558 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 559 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 560 PetscFunctionReturn(0); 561 } 562 563 #undef __FUNCT__ 564 #define __FUNCT__ "MatView_SeqAIJ" 565 int MatView_SeqAIJ(Mat A,PetscViewer viewer) 566 { 567 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 568 int ierr; 569 PetscTruth issocket,isascii,isbinary,isdraw; 570 571 PetscFunctionBegin; 572 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 573 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 574 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 575 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 576 if (issocket) { 577 ierr = PetscViewerSocketPutSparse_Private(viewer,A->m,A->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 578 } else if (isascii) { 579 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 580 } else if (isbinary) { 581 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 582 } else if (isdraw) { 583 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 584 } else { 585 SETERRQ1(1,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 586 } 587 PetscFunctionReturn(0); 588 } 589 590 EXTERN int Mat_AIJ_CheckInode(Mat,PetscTruth); 591 #undef __FUNCT__ 592 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 593 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 594 { 595 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 596 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax,ierr; 597 int m = A->m,*ip,N,*ailen = a->ilen,rmax = 0; 598 PetscScalar *aa = a->a,*ap; 599 600 PetscFunctionBegin; 601 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 602 603 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 604 for (i=1; i<m; i++) { 605 /* move each row back by the amount of empty slots (fshift) before it*/ 606 fshift += imax[i-1] - ailen[i-1]; 607 rmax = PetscMax(rmax,ailen[i]); 608 if (fshift) { 609 ip = aj + ai[i] ; 610 ap = aa + ai[i] ; 611 N = ailen[i]; 612 for (j=0; j<N; j++) { 613 ip[j-fshift] = ip[j]; 614 ap[j-fshift] = ap[j]; 615 } 616 } 617 ai[i] = ai[i-1] + ailen[i-1]; 618 } 619 if (m) { 620 fshift += imax[m-1] - ailen[m-1]; 621 ai[m] = ai[m-1] + ailen[m-1]; 622 } 623 /* reset ilen and imax for each row */ 624 for (i=0; i<m; i++) { 625 ailen[i] = imax[i] = ai[i+1] - ai[i]; 626 } 627 a->nz = ai[m]; 628 629 /* diagonals may have moved, so kill the diagonal pointers */ 630 if (fshift && a->diag) { 631 ierr = PetscFree(a->diag);CHKERRQ(ierr); 632 PetscLogObjectMemory(A,-(m+1)*sizeof(int)); 633 a->diag = 0; 634 } 635 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded,%d used\n",m,A->n,fshift,a->nz); 636 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %d\n",a->reallocs); 637 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 638 a->reallocs = 0; 639 A->info.nz_unneeded = (double)fshift; 640 a->rmax = rmax; 641 642 /* check out for identical nodes. If found, use inode functions */ 643 ierr = Mat_AIJ_CheckInode(A,(PetscTruth)(!fshift));CHKERRQ(ierr); 644 645 PetscFunctionReturn(0); 646 } 647 648 #undef __FUNCT__ 649 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 650 int MatZeroEntries_SeqAIJ(Mat A) 651 { 652 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 653 int ierr; 654 655 PetscFunctionBegin; 656 ierr = PetscMemzero(a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 657 PetscFunctionReturn(0); 658 } 659 660 #undef __FUNCT__ 661 #define __FUNCT__ "MatDestroy_SeqAIJ" 662 int MatDestroy_SeqAIJ(Mat A) 663 { 664 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 665 int ierr; 666 667 PetscFunctionBegin; 668 #if defined(PETSC_USE_LOG) 669 PetscLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",A->m,A->n,a->nz); 670 #endif 671 if (a->freedata) { 672 ierr = PetscFree(a->a);CHKERRQ(ierr); 673 if (!a->singlemalloc) { 674 ierr = PetscFree(a->i);CHKERRQ(ierr); 675 ierr = PetscFree(a->j);CHKERRQ(ierr); 676 } 677 } 678 if (a->row) { 679 ierr = ISDestroy(a->row);CHKERRQ(ierr); 680 } 681 if (a->col) { 682 ierr = ISDestroy(a->col);CHKERRQ(ierr); 683 } 684 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 685 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 686 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 687 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 688 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 689 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 690 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 691 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 692 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 693 if (a->xtoy) {ierr = PetscFree(a->xtoy);CHKERRQ(ierr);} 694 695 ierr = PetscFree(a);CHKERRQ(ierr); 696 PetscFunctionReturn(0); 697 } 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "MatCompress_SeqAIJ" 701 int MatCompress_SeqAIJ(Mat A) 702 { 703 PetscFunctionBegin; 704 PetscFunctionReturn(0); 705 } 706 707 #undef __FUNCT__ 708 #define __FUNCT__ "MatSetOption_SeqAIJ" 709 int MatSetOption_SeqAIJ(Mat A,MatOption op) 710 { 711 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 712 713 PetscFunctionBegin; 714 switch (op) { 715 case MAT_ROW_ORIENTED: 716 a->roworiented = PETSC_TRUE; 717 break; 718 case MAT_KEEP_ZEROED_ROWS: 719 a->keepzeroedrows = PETSC_TRUE; 720 break; 721 case MAT_COLUMN_ORIENTED: 722 a->roworiented = PETSC_FALSE; 723 break; 724 case MAT_COLUMNS_SORTED: 725 a->sorted = PETSC_TRUE; 726 break; 727 case MAT_COLUMNS_UNSORTED: 728 a->sorted = PETSC_FALSE; 729 break; 730 case MAT_NO_NEW_NONZERO_LOCATIONS: 731 a->nonew = 1; 732 break; 733 case MAT_NEW_NONZERO_LOCATION_ERR: 734 a->nonew = -1; 735 break; 736 case MAT_NEW_NONZERO_ALLOCATION_ERR: 737 a->nonew = -2; 738 break; 739 case MAT_YES_NEW_NONZERO_LOCATIONS: 740 a->nonew = 0; 741 break; 742 case MAT_IGNORE_ZERO_ENTRIES: 743 a->ignorezeroentries = PETSC_TRUE; 744 break; 745 case MAT_USE_INODES: 746 a->inode.use = PETSC_TRUE; 747 break; 748 case MAT_DO_NOT_USE_INODES: 749 a->inode.use = PETSC_FALSE; 750 break; 751 case MAT_ROWS_SORTED: 752 case MAT_ROWS_UNSORTED: 753 case MAT_YES_NEW_DIAGONALS: 754 case MAT_IGNORE_OFF_PROC_ENTRIES: 755 case MAT_USE_HASH_TABLE: 756 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 757 break; 758 case MAT_NO_NEW_DIAGONALS: 759 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 760 case MAT_INODE_LIMIT_1: 761 a->inode.limit = 1; 762 break; 763 case MAT_INODE_LIMIT_2: 764 a->inode.limit = 2; 765 break; 766 case MAT_INODE_LIMIT_3: 767 a->inode.limit = 3; 768 break; 769 case MAT_INODE_LIMIT_4: 770 a->inode.limit = 4; 771 break; 772 case MAT_INODE_LIMIT_5: 773 a->inode.limit = 5; 774 break; 775 default: 776 SETERRQ(PETSC_ERR_SUP,"unknown option"); 777 } 778 PetscFunctionReturn(0); 779 } 780 781 #undef __FUNCT__ 782 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 783 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 784 { 785 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 786 int i,j,n,ierr; 787 PetscScalar *x,zero = 0.0; 788 789 PetscFunctionBegin; 790 ierr = VecSet(&zero,v);CHKERRQ(ierr); 791 ierr = VecGetArrayFast(v,&x);CHKERRQ(ierr); 792 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 793 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 794 for (i=0; i<A->m; i++) { 795 for (j=a->i[i]; j<a->i[i+1]; j++) { 796 if (a->j[j] == i) { 797 x[i] = a->a[j]; 798 break; 799 } 800 } 801 } 802 ierr = VecRestoreArrayFast(v,&x);CHKERRQ(ierr); 803 PetscFunctionReturn(0); 804 } 805 806 807 #undef __FUNCT__ 808 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 809 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 810 { 811 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 812 PetscScalar *x,*y; 813 int ierr,m = A->m; 814 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 815 PetscScalar *v,alpha; 816 int n,i,*idx; 817 #endif 818 819 PetscFunctionBegin; 820 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 821 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 822 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 823 824 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 825 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 826 #else 827 for (i=0; i<m; i++) { 828 idx = a->j + a->i[i] ; 829 v = a->a + a->i[i] ; 830 n = a->i[i+1] - a->i[i]; 831 alpha = x[i]; 832 while (n-->0) {y[*idx++] += alpha * *v++;} 833 } 834 #endif 835 PetscLogFlops(2*a->nz); 836 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 837 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 838 PetscFunctionReturn(0); 839 } 840 841 #undef __FUNCT__ 842 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 843 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 844 { 845 PetscScalar zero = 0.0; 846 int ierr; 847 848 PetscFunctionBegin; 849 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 850 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 851 PetscFunctionReturn(0); 852 } 853 854 855 #undef __FUNCT__ 856 #define __FUNCT__ "MatMult_SeqAIJ" 857 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 858 { 859 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 860 PetscScalar *x,*y,*v; 861 int ierr,m = A->m,*idx,*ii; 862 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 863 int n,i,jrow,j; 864 PetscScalar sum; 865 #endif 866 867 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 868 #pragma disjoint(*x,*y,*v) 869 #endif 870 871 PetscFunctionBegin; 872 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 873 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 874 idx = a->j; 875 v = a->a; 876 ii = a->i; 877 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 878 fortranmultaij_(&m,x,ii,idx,v,y); 879 #else 880 for (i=0; i<m; i++) { 881 jrow = ii[i]; 882 n = ii[i+1] - jrow; 883 sum = 0.0; 884 for (j=0; j<n; j++) { 885 sum += v[jrow]*x[idx[jrow]]; jrow++; 886 } 887 y[i] = sum; 888 } 889 #endif 890 PetscLogFlops(2*a->nz - m); 891 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 892 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 893 PetscFunctionReturn(0); 894 } 895 896 #undef __FUNCT__ 897 #define __FUNCT__ "MatMultAdd_SeqAIJ" 898 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 899 { 900 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 901 PetscScalar *x,*y,*z,*v; 902 int ierr,m = A->m,*idx,*ii; 903 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 904 int n,i,jrow,j; 905 PetscScalar sum; 906 #endif 907 908 PetscFunctionBegin; 909 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 910 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 911 if (zz != yy) { 912 ierr = VecGetArrayFast(zz,&z);CHKERRQ(ierr); 913 } else { 914 z = y; 915 } 916 917 idx = a->j; 918 v = a->a; 919 ii = a->i; 920 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 921 fortranmultaddaij_(&m,x,ii,idx,v,y,z); 922 #else 923 for (i=0; i<m; i++) { 924 jrow = ii[i]; 925 n = ii[i+1] - jrow; 926 sum = y[i]; 927 for (j=0; j<n; j++) { 928 sum += v[jrow]*x[idx[jrow]]; jrow++; 929 } 930 z[i] = sum; 931 } 932 #endif 933 PetscLogFlops(2*a->nz); 934 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 935 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 936 if (zz != yy) { 937 ierr = VecRestoreArrayFast(zz,&z);CHKERRQ(ierr); 938 } 939 PetscFunctionReturn(0); 940 } 941 942 /* 943 Adds diagonal pointers to sparse matrix structure. 944 */ 945 #undef __FUNCT__ 946 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 947 int MatMarkDiagonal_SeqAIJ(Mat A) 948 { 949 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 950 int i,j,*diag,m = A->m,ierr; 951 952 PetscFunctionBegin; 953 if (a->diag) PetscFunctionReturn(0); 954 955 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 956 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 957 for (i=0; i<A->m; i++) { 958 diag[i] = a->i[i+1]; 959 for (j=a->i[i]; j<a->i[i+1]; j++) { 960 if (a->j[j] == i) { 961 diag[i] = j; 962 break; 963 } 964 } 965 } 966 a->diag = diag; 967 PetscFunctionReturn(0); 968 } 969 970 /* 971 Checks for missing diagonals 972 */ 973 #undef __FUNCT__ 974 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 975 int MatMissingDiagonal_SeqAIJ(Mat A) 976 { 977 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 978 int *diag,*jj = a->j,i,ierr; 979 980 PetscFunctionBegin; 981 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 982 diag = a->diag; 983 for (i=0; i<A->m; i++) { 984 if (jj[diag[i]] != i) { 985 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 986 } 987 } 988 PetscFunctionReturn(0); 989 } 990 991 #undef __FUNCT__ 992 #define __FUNCT__ "MatRelax_SeqAIJ" 993 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx) 994 { 995 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 996 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 997 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 998 int ierr,n = A->n,m = A->m,i; 999 const int *idx,*diag; 1000 1001 PetscFunctionBegin; 1002 its = its*lits; 1003 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits); 1004 1005 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1006 diag = a->diag; 1007 if (!a->idiag) { 1008 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1009 a->ssor = a->idiag + m; 1010 mdiag = a->ssor + m; 1011 1012 v = a->a; 1013 1014 /* this is wrong when fshift omega changes each iteration */ 1015 if (omega == 1.0 && fshift == 0.0) { 1016 for (i=0; i<m; i++) { 1017 mdiag[i] = v[diag[i]]; 1018 a->idiag[i] = 1.0/v[diag[i]]; 1019 } 1020 PetscLogFlops(m); 1021 } else { 1022 for (i=0; i<m; i++) { 1023 mdiag[i] = v[diag[i]]; 1024 a->idiag[i] = omega/(fshift + v[diag[i]]); 1025 } 1026 PetscLogFlops(2*m); 1027 } 1028 } 1029 t = a->ssor; 1030 idiag = a->idiag; 1031 mdiag = a->idiag + 2*m; 1032 1033 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 1034 if (xx != bb) { 1035 ierr = VecGetArrayFast(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1036 } else { 1037 b = x; 1038 } 1039 1040 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1041 xs = x; 1042 if (flag == SOR_APPLY_UPPER) { 1043 /* apply (U + D/omega) to the vector */ 1044 bs = b; 1045 for (i=0; i<m; i++) { 1046 d = fshift + a->a[diag[i]]; 1047 n = a->i[i+1] - diag[i] - 1; 1048 idx = a->j + diag[i] + 1; 1049 v = a->a + diag[i] + 1; 1050 sum = b[i]*d/omega; 1051 SPARSEDENSEDOT(sum,bs,v,idx,n); 1052 x[i] = sum; 1053 } 1054 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1055 if (bb != xx) {ierr = VecRestoreArrayFast(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1056 PetscLogFlops(a->nz); 1057 PetscFunctionReturn(0); 1058 } 1059 1060 1061 /* Let A = L + U + D; where L is lower trianglar, 1062 U is upper triangular, E is diagonal; This routine applies 1063 1064 (L + E)^{-1} A (U + E)^{-1} 1065 1066 to a vector efficiently using Eisenstat's trick. This is for 1067 the case of SSOR preconditioner, so E is D/omega where omega 1068 is the relaxation factor. 1069 */ 1070 1071 if (flag == SOR_APPLY_LOWER) { 1072 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1073 } else if (flag & SOR_EISENSTAT) { 1074 /* Let A = L + U + D; where L is lower trianglar, 1075 U is upper triangular, E is diagonal; This routine applies 1076 1077 (L + E)^{-1} A (U + E)^{-1} 1078 1079 to a vector efficiently using Eisenstat's trick. This is for 1080 the case of SSOR preconditioner, so E is D/omega where omega 1081 is the relaxation factor. 1082 */ 1083 scale = (2.0/omega) - 1.0; 1084 1085 /* x = (E + U)^{-1} b */ 1086 for (i=m-1; i>=0; i--) { 1087 n = a->i[i+1] - diag[i] - 1; 1088 idx = a->j + diag[i] + 1; 1089 v = a->a + diag[i] + 1; 1090 sum = b[i]; 1091 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1092 x[i] = sum*idiag[i]; 1093 } 1094 1095 /* t = b - (2*E - D)x */ 1096 v = a->a; 1097 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1098 1099 /* t = (E + L)^{-1}t */ 1100 ts = t; 1101 diag = a->diag; 1102 for (i=0; i<m; i++) { 1103 n = diag[i] - a->i[i]; 1104 idx = a->j + a->i[i]; 1105 v = a->a + a->i[i]; 1106 sum = t[i]; 1107 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1108 t[i] = sum*idiag[i]; 1109 /* x = x + t */ 1110 x[i] += t[i]; 1111 } 1112 1113 PetscLogFlops(6*m-1 + 2*a->nz); 1114 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1115 if (bb != xx) {ierr = VecRestoreArrayFast(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1116 PetscFunctionReturn(0); 1117 } 1118 if (flag & SOR_ZERO_INITIAL_GUESS) { 1119 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1120 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1121 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(int *)diag,idiag,a->a,(void*)b); 1122 #else 1123 for (i=0; i<m; i++) { 1124 n = diag[i] - a->i[i]; 1125 idx = a->j + a->i[i]; 1126 v = a->a + a->i[i]; 1127 sum = b[i]; 1128 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1129 x[i] = sum*idiag[i]; 1130 } 1131 #endif 1132 xb = x; 1133 PetscLogFlops(a->nz); 1134 } else xb = b; 1135 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1136 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1137 for (i=0; i<m; i++) { 1138 x[i] *= mdiag[i]; 1139 } 1140 PetscLogFlops(m); 1141 } 1142 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1143 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1144 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(int*)diag,idiag,a->a,(void*)xb); 1145 #else 1146 for (i=m-1; i>=0; i--) { 1147 n = a->i[i+1] - diag[i] - 1; 1148 idx = a->j + diag[i] + 1; 1149 v = a->a + diag[i] + 1; 1150 sum = xb[i]; 1151 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1152 x[i] = sum*idiag[i]; 1153 } 1154 #endif 1155 PetscLogFlops(a->nz); 1156 } 1157 its--; 1158 } 1159 while (its--) { 1160 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1161 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1162 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(int*)diag,a->a,(void*)b); 1163 #else 1164 for (i=0; i<m; i++) { 1165 d = fshift + a->a[diag[i]]; 1166 n = a->i[i+1] - a->i[i]; 1167 idx = a->j + a->i[i]; 1168 v = a->a + a->i[i]; 1169 sum = b[i]; 1170 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1171 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1172 } 1173 #endif 1174 PetscLogFlops(a->nz); 1175 } 1176 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1177 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1178 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(int*)diag,a->a,(void*)b); 1179 #else 1180 for (i=m-1; i>=0; i--) { 1181 d = fshift + a->a[diag[i]]; 1182 n = a->i[i+1] - a->i[i]; 1183 idx = a->j + a->i[i]; 1184 v = a->a + a->i[i]; 1185 sum = b[i]; 1186 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1187 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1188 } 1189 #endif 1190 PetscLogFlops(a->nz); 1191 } 1192 } 1193 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1194 if (bb != xx) {ierr = VecRestoreArrayFast(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNCT__ 1199 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1200 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1201 { 1202 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1203 1204 PetscFunctionBegin; 1205 info->rows_global = (double)A->m; 1206 info->columns_global = (double)A->n; 1207 info->rows_local = (double)A->m; 1208 info->columns_local = (double)A->n; 1209 info->block_size = 1.0; 1210 info->nz_allocated = (double)a->maxnz; 1211 info->nz_used = (double)a->nz; 1212 info->nz_unneeded = (double)(a->maxnz - a->nz); 1213 info->assemblies = (double)A->num_ass; 1214 info->mallocs = (double)a->reallocs; 1215 info->memory = A->mem; 1216 if (A->factor) { 1217 info->fill_ratio_given = A->info.fill_ratio_given; 1218 info->fill_ratio_needed = A->info.fill_ratio_needed; 1219 info->factor_mallocs = A->info.factor_mallocs; 1220 } else { 1221 info->fill_ratio_given = 0; 1222 info->fill_ratio_needed = 0; 1223 info->factor_mallocs = 0; 1224 } 1225 PetscFunctionReturn(0); 1226 } 1227 1228 EXTERN int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,MatFactorInfo*,Mat*); 1229 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1230 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatFactorInfo*); 1231 EXTERN int MatSolve_SeqAIJ(Mat,Vec,Vec); 1232 EXTERN int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1233 EXTERN int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1234 EXTERN int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1235 1236 #undef __FUNCT__ 1237 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1238 int MatZeroRows_SeqAIJ(Mat A,IS is,const PetscScalar *diag) 1239 { 1240 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1241 int i,ierr,N,*rows,m = A->m - 1; 1242 1243 PetscFunctionBegin; 1244 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1245 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1246 if (a->keepzeroedrows) { 1247 for (i=0; i<N; i++) { 1248 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1249 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1250 } 1251 if (diag) { 1252 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1253 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1254 for (i=0; i<N; i++) { 1255 a->a[a->diag[rows[i]]] = *diag; 1256 } 1257 } 1258 } else { 1259 if (diag) { 1260 for (i=0; i<N; i++) { 1261 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1262 if (a->ilen[rows[i]] > 0) { 1263 a->ilen[rows[i]] = 1; 1264 a->a[a->i[rows[i]]] = *diag; 1265 a->j[a->i[rows[i]]] = rows[i]; 1266 } else { /* in case row was completely empty */ 1267 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1268 } 1269 } 1270 } else { 1271 for (i=0; i<N; i++) { 1272 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1273 a->ilen[rows[i]] = 0; 1274 } 1275 } 1276 } 1277 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1278 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1279 PetscFunctionReturn(0); 1280 } 1281 1282 #undef __FUNCT__ 1283 #define __FUNCT__ "MatGetRow_SeqAIJ" 1284 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1285 { 1286 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1287 int *itmp; 1288 1289 PetscFunctionBegin; 1290 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1291 1292 *nz = a->i[row+1] - a->i[row]; 1293 if (v) *v = a->a + a->i[row]; 1294 if (idx) { 1295 itmp = a->j + a->i[row]; 1296 if (*nz) { 1297 *idx = itmp; 1298 } 1299 else *idx = 0; 1300 } 1301 PetscFunctionReturn(0); 1302 } 1303 1304 /* remove this function? */ 1305 #undef __FUNCT__ 1306 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1307 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1308 { 1309 /* Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1310 int ierr; */ 1311 1312 PetscFunctionBegin; 1313 /* if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} */ 1314 PetscFunctionReturn(0); 1315 } 1316 1317 #undef __FUNCT__ 1318 #define __FUNCT__ "MatNorm_SeqAIJ" 1319 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1320 { 1321 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1322 PetscScalar *v = a->a; 1323 PetscReal sum = 0.0; 1324 int i,j,ierr; 1325 1326 PetscFunctionBegin; 1327 if (type == NORM_FROBENIUS) { 1328 for (i=0; i<a->nz; i++) { 1329 #if defined(PETSC_USE_COMPLEX) 1330 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1331 #else 1332 sum += (*v)*(*v); v++; 1333 #endif 1334 } 1335 *nrm = sqrt(sum); 1336 } else if (type == NORM_1) { 1337 PetscReal *tmp; 1338 int *jj = a->j; 1339 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1340 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1341 *nrm = 0.0; 1342 for (j=0; j<a->nz; j++) { 1343 tmp[*jj++] += PetscAbsScalar(*v); v++; 1344 } 1345 for (j=0; j<A->n; j++) { 1346 if (tmp[j] > *nrm) *nrm = tmp[j]; 1347 } 1348 ierr = PetscFree(tmp);CHKERRQ(ierr); 1349 } else if (type == NORM_INFINITY) { 1350 *nrm = 0.0; 1351 for (j=0; j<A->m; j++) { 1352 v = a->a + a->i[j]; 1353 sum = 0.0; 1354 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1355 sum += PetscAbsScalar(*v); v++; 1356 } 1357 if (sum > *nrm) *nrm = sum; 1358 } 1359 } else { 1360 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1361 } 1362 PetscFunctionReturn(0); 1363 } 1364 1365 #undef __FUNCT__ 1366 #define __FUNCT__ "MatTranspose_SeqAIJ" 1367 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1368 { 1369 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1370 Mat C; 1371 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1372 PetscScalar *array = a->a; 1373 1374 PetscFunctionBegin; 1375 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1376 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1377 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1378 1379 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1380 ierr = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1381 ierr = PetscFree(col);CHKERRQ(ierr); 1382 for (i=0; i<m; i++) { 1383 len = ai[i+1]-ai[i]; 1384 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1385 array += len; 1386 aj += len; 1387 } 1388 1389 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1390 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1391 1392 if (B) { 1393 *B = C; 1394 } else { 1395 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1396 } 1397 PetscFunctionReturn(0); 1398 } 1399 1400 EXTERN_C_BEGIN 1401 #undef __FUNCT__ 1402 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1403 int MatIsSymmetric_SeqAIJ(Mat A,Mat B,PetscTruth *f) 1404 { 1405 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1406 int *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1407 int ma,na,mb,nb, i,ierr; 1408 1409 PetscFunctionBegin; 1410 bij = (Mat_SeqAIJ *) B->data; 1411 1412 ierr = MatGetSize(A,&ma,&na); CHKERRQ(ierr); 1413 ierr = MatGetSize(B,&mb,&nb); CHKERRQ(ierr); 1414 if (ma!=nb || na!=mb) 1415 SETERRQ(1,"Incompatible A/B sizes for symmetry test"); 1416 aii = aij->i; bii = bij->i; 1417 adx = aij->j; bdx = bij->j; 1418 va = aij->a; vb = bij->a; 1419 ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr); 1420 ierr = PetscMalloc(mb*sizeof(int),&bptr); CHKERRQ(ierr); 1421 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1422 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1423 1424 *f = PETSC_TRUE; 1425 for (i=0; i<ma; i++) { 1426 /*printf("row %d spans %d--%d; we start @ %d\n", 1427 i,idx[ii[i]],idx[ii[i+1]-1],idx[aptr[i]]);*/ 1428 while (aptr[i]<aii[i+1]) { 1429 int idc,idr; PetscScalar vc,vr; 1430 /* column/row index/value */ 1431 idc = adx[aptr[i]]; idr = bdx[bptr[idc]]; 1432 vc = va[aptr[i]]; vr = vb[bptr[idc]]; 1433 /*printf("comparing %d: (%d,%d)=%e to (%d,%d)=%e\n", 1434 aptr[i],i,idc,vc,idc,idr,vr);*/ 1435 if (i!=idr || vc!=vr) { 1436 *f = PETSC_FALSE; goto done; 1437 } else { 1438 aptr[i]++; if (B || i!=idc) bptr[idc]++; 1439 } 1440 } 1441 } 1442 done: 1443 ierr = PetscFree(aptr); CHKERRQ(ierr); 1444 if (B) { 1445 ierr = PetscFree(bptr); CHKERRQ(ierr); 1446 } 1447 1448 PetscFunctionReturn(0); 1449 } 1450 EXTERN_C_END 1451 1452 #undef __FUNCT__ 1453 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1454 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1455 { 1456 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1457 PetscScalar *l,*r,x,*v; 1458 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1459 1460 PetscFunctionBegin; 1461 if (ll) { 1462 /* The local size is used so that VecMPI can be passed to this routine 1463 by MatDiagonalScale_MPIAIJ */ 1464 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1465 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1466 ierr = VecGetArrayFast(ll,&l);CHKERRQ(ierr); 1467 v = a->a; 1468 for (i=0; i<m; i++) { 1469 x = l[i]; 1470 M = a->i[i+1] - a->i[i]; 1471 for (j=0; j<M; j++) { (*v++) *= x;} 1472 } 1473 ierr = VecRestoreArrayFast(ll,&l);CHKERRQ(ierr); 1474 PetscLogFlops(nz); 1475 } 1476 if (rr) { 1477 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1478 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1479 ierr = VecGetArrayFast(rr,&r);CHKERRQ(ierr); 1480 v = a->a; jj = a->j; 1481 for (i=0; i<nz; i++) { 1482 (*v++) *= r[*jj++]; 1483 } 1484 ierr = VecRestoreArrayFast(rr,&r);CHKERRQ(ierr); 1485 PetscLogFlops(nz); 1486 } 1487 PetscFunctionReturn(0); 1488 } 1489 1490 #undef __FUNCT__ 1491 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1492 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1493 { 1494 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1495 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1496 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1497 int *irow,*icol,nrows,ncols; 1498 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1499 PetscScalar *a_new,*mat_a; 1500 Mat C; 1501 PetscTruth stride; 1502 1503 PetscFunctionBegin; 1504 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1505 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1506 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1507 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1508 1509 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1510 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1511 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1512 1513 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1514 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1515 if (stride && step == 1) { 1516 /* special case of contiguous rows */ 1517 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1518 starts = lens + nrows; 1519 /* loop over new rows determining lens and starting points */ 1520 for (i=0; i<nrows; i++) { 1521 kstart = ai[irow[i]]; 1522 kend = kstart + ailen[irow[i]]; 1523 for (k=kstart; k<kend; k++) { 1524 if (aj[k] >= first) { 1525 starts[i] = k; 1526 break; 1527 } 1528 } 1529 sum = 0; 1530 while (k < kend) { 1531 if (aj[k++] >= first+ncols) break; 1532 sum++; 1533 } 1534 lens[i] = sum; 1535 } 1536 /* create submatrix */ 1537 if (scall == MAT_REUSE_MATRIX) { 1538 int n_cols,n_rows; 1539 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1540 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1541 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1542 C = *B; 1543 } else { 1544 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1545 } 1546 c = (Mat_SeqAIJ*)C->data; 1547 1548 /* loop over rows inserting into submatrix */ 1549 a_new = c->a; 1550 j_new = c->j; 1551 i_new = c->i; 1552 1553 for (i=0; i<nrows; i++) { 1554 ii = starts[i]; 1555 lensi = lens[i]; 1556 for (k=0; k<lensi; k++) { 1557 *j_new++ = aj[ii+k] - first; 1558 } 1559 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1560 a_new += lensi; 1561 i_new[i+1] = i_new[i] + lensi; 1562 c->ilen[i] = lensi; 1563 } 1564 ierr = PetscFree(lens);CHKERRQ(ierr); 1565 } else { 1566 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1567 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1568 1569 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1570 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1571 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1572 /* determine lens of each row */ 1573 for (i=0; i<nrows; i++) { 1574 kstart = ai[irow[i]]; 1575 kend = kstart + a->ilen[irow[i]]; 1576 lens[i] = 0; 1577 for (k=kstart; k<kend; k++) { 1578 if (smap[aj[k]]) { 1579 lens[i]++; 1580 } 1581 } 1582 } 1583 /* Create and fill new matrix */ 1584 if (scall == MAT_REUSE_MATRIX) { 1585 PetscTruth equal; 1586 1587 c = (Mat_SeqAIJ *)((*B)->data); 1588 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1589 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1590 if (!equal) { 1591 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1592 } 1593 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1594 C = *B; 1595 } else { 1596 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1597 } 1598 c = (Mat_SeqAIJ *)(C->data); 1599 for (i=0; i<nrows; i++) { 1600 row = irow[i]; 1601 kstart = ai[row]; 1602 kend = kstart + a->ilen[row]; 1603 mat_i = c->i[i]; 1604 mat_j = c->j + mat_i; 1605 mat_a = c->a + mat_i; 1606 mat_ilen = c->ilen + i; 1607 for (k=kstart; k<kend; k++) { 1608 if ((tcol=smap[a->j[k]])) { 1609 *mat_j++ = tcol - 1; 1610 *mat_a++ = a->a[k]; 1611 (*mat_ilen)++; 1612 1613 } 1614 } 1615 } 1616 /* Free work space */ 1617 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1618 ierr = PetscFree(smap);CHKERRQ(ierr); 1619 ierr = PetscFree(lens);CHKERRQ(ierr); 1620 } 1621 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1622 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1623 1624 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1625 *B = C; 1626 PetscFunctionReturn(0); 1627 } 1628 1629 /* 1630 */ 1631 #undef __FUNCT__ 1632 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1633 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1634 { 1635 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1636 int ierr; 1637 Mat outA; 1638 PetscTruth row_identity,col_identity; 1639 1640 PetscFunctionBegin; 1641 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1642 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1643 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1644 if (!row_identity || !col_identity) { 1645 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1646 } 1647 1648 outA = inA; 1649 inA->factor = FACTOR_LU; 1650 a->row = row; 1651 a->col = col; 1652 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1653 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1654 1655 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1656 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1657 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1658 PetscLogObjectParent(inA,a->icol); 1659 1660 if (!a->solve_work) { /* this matrix may have been factored before */ 1661 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1662 } 1663 1664 if (!a->diag) { 1665 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1666 } 1667 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1668 PetscFunctionReturn(0); 1669 } 1670 1671 #include "petscblaslapack.h" 1672 #undef __FUNCT__ 1673 #define __FUNCT__ "MatScale_SeqAIJ" 1674 int MatScale_SeqAIJ(const PetscScalar *alpha,Mat inA) 1675 { 1676 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1677 int one = 1; 1678 1679 PetscFunctionBegin; 1680 BLscal_(&a->nz,(PetscScalar*)alpha,a->a,&one); 1681 PetscLogFlops(a->nz); 1682 PetscFunctionReturn(0); 1683 } 1684 1685 #undef __FUNCT__ 1686 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1687 int MatGetSubMatrices_SeqAIJ(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1688 { 1689 int ierr,i; 1690 1691 PetscFunctionBegin; 1692 if (scall == MAT_INITIAL_MATRIX) { 1693 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1694 } 1695 1696 for (i=0; i<n; i++) { 1697 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1698 } 1699 PetscFunctionReturn(0); 1700 } 1701 1702 #undef __FUNCT__ 1703 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1704 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1705 { 1706 PetscFunctionBegin; 1707 *bs = 1; 1708 PetscFunctionReturn(0); 1709 } 1710 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1713 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS is[],int ov) 1714 { 1715 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1716 int row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1717 int start,end,*ai,*aj; 1718 PetscBT table; 1719 1720 PetscFunctionBegin; 1721 m = A->m; 1722 ai = a->i; 1723 aj = a->j; 1724 1725 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1726 1727 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1728 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1729 1730 for (i=0; i<is_max; i++) { 1731 /* Initialize the two local arrays */ 1732 isz = 0; 1733 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1734 1735 /* Extract the indices, assume there can be duplicate entries */ 1736 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1737 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1738 1739 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1740 for (j=0; j<n ; ++j){ 1741 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1742 } 1743 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1744 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1745 1746 k = 0; 1747 for (j=0; j<ov; j++){ /* for each overlap */ 1748 n = isz; 1749 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1750 row = nidx[k]; 1751 start = ai[row]; 1752 end = ai[row+1]; 1753 for (l = start; l<end ; l++){ 1754 val = aj[l] ; 1755 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1756 } 1757 } 1758 } 1759 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1760 } 1761 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1762 ierr = PetscFree(nidx);CHKERRQ(ierr); 1763 PetscFunctionReturn(0); 1764 } 1765 1766 /* -------------------------------------------------------------- */ 1767 #undef __FUNCT__ 1768 #define __FUNCT__ "MatPermute_SeqAIJ" 1769 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1770 { 1771 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1772 PetscScalar *vwork; 1773 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1774 int *row,*col,*cnew,j,*lens; 1775 IS icolp,irowp; 1776 1777 PetscFunctionBegin; 1778 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1779 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1780 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1781 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1782 1783 /* determine lengths of permuted rows */ 1784 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1785 for (i=0; i<m; i++) { 1786 lens[row[i]] = a->i[i+1] - a->i[i]; 1787 } 1788 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1789 ierr = PetscFree(lens);CHKERRQ(ierr); 1790 1791 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1792 for (i=0; i<m; i++) { 1793 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1794 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1795 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1796 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1797 } 1798 ierr = PetscFree(cnew);CHKERRQ(ierr); 1799 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1800 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1801 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1802 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1803 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1804 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1805 PetscFunctionReturn(0); 1806 } 1807 1808 #undef __FUNCT__ 1809 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1810 int MatPrintHelp_SeqAIJ(Mat A) 1811 { 1812 static PetscTruth called = PETSC_FALSE; 1813 MPI_Comm comm = A->comm; 1814 int ierr; 1815 1816 PetscFunctionBegin; 1817 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1818 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1819 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1820 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1821 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1822 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1823 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1824 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1825 #endif 1826 PetscFunctionReturn(0); 1827 } 1828 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1829 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1830 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatFactorInfo*,IS,IS,Mat*); 1831 #undef __FUNCT__ 1832 #define __FUNCT__ "MatCopy_SeqAIJ" 1833 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1834 { 1835 int ierr; 1836 1837 PetscFunctionBegin; 1838 /* If the two matrices have the same copy implementation, use fast copy. */ 1839 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1840 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1841 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1842 1843 if (a->i[A->m] != b->i[B->m]) { 1844 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1845 } 1846 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1847 } else { 1848 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1849 } 1850 PetscFunctionReturn(0); 1851 } 1852 1853 #undef __FUNCT__ 1854 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1855 int MatSetUpPreallocation_SeqAIJ(Mat A) 1856 { 1857 int ierr; 1858 1859 PetscFunctionBegin; 1860 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1861 PetscFunctionReturn(0); 1862 } 1863 1864 #undef __FUNCT__ 1865 #define __FUNCT__ "MatGetArray_SeqAIJ" 1866 int MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1867 { 1868 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1869 PetscFunctionBegin; 1870 *array = a->a; 1871 PetscFunctionReturn(0); 1872 } 1873 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1876 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1877 { 1878 PetscFunctionBegin; 1879 PetscFunctionReturn(0); 1880 } 1881 1882 #undef __FUNCT__ 1883 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1884 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1885 { 1886 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1887 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1888 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1889 PetscScalar *vscale_array; 1890 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1891 Vec w1,w2,w3; 1892 void *fctx = coloring->fctx; 1893 PetscTruth flg; 1894 1895 PetscFunctionBegin; 1896 if (!coloring->w1) { 1897 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1898 PetscLogObjectParent(coloring,coloring->w1); 1899 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1900 PetscLogObjectParent(coloring,coloring->w2); 1901 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1902 PetscLogObjectParent(coloring,coloring->w3); 1903 } 1904 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1905 1906 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1907 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1908 if (flg) { 1909 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1910 } else { 1911 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1912 } 1913 1914 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1915 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1916 1917 /* 1918 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1919 coloring->F for the coarser grids from the finest 1920 */ 1921 if (coloring->F) { 1922 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1923 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1924 if (m1 != m2) { 1925 coloring->F = 0; 1926 } 1927 } 1928 1929 if (coloring->F) { 1930 w1 = coloring->F; 1931 coloring->F = 0; 1932 } else { 1933 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1934 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1935 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1936 } 1937 1938 /* 1939 Compute all the scale factors and share with other processors 1940 */ 1941 ierr = VecGetArrayFast(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1942 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1943 for (k=0; k<coloring->ncolors; k++) { 1944 /* 1945 Loop over each column associated with color adding the 1946 perturbation to the vector w3. 1947 */ 1948 for (l=0; l<coloring->ncolumns[k]; l++) { 1949 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1950 dx = xx[col]; 1951 if (dx == 0.0) dx = 1.0; 1952 #if !defined(PETSC_USE_COMPLEX) 1953 if (dx < umin && dx >= 0.0) dx = umin; 1954 else if (dx < 0.0 && dx > -umin) dx = -umin; 1955 #else 1956 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1957 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1958 #endif 1959 dx *= epsilon; 1960 vscale_array[col] = 1.0/dx; 1961 } 1962 } 1963 vscale_array = vscale_array + start;ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1964 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1965 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1966 1967 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1968 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1969 1970 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1971 else vscaleforrow = coloring->columnsforrow; 1972 1973 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1974 /* 1975 Loop over each color 1976 */ 1977 for (k=0; k<coloring->ncolors; k++) { 1978 coloring->currentcolor = k; 1979 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 1980 ierr = VecGetArrayFast(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 1981 /* 1982 Loop over each column associated with color adding the 1983 perturbation to the vector w3. 1984 */ 1985 for (l=0; l<coloring->ncolumns[k]; l++) { 1986 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1987 dx = xx[col]; 1988 if (dx == 0.0) dx = 1.0; 1989 #if !defined(PETSC_USE_COMPLEX) 1990 if (dx < umin && dx >= 0.0) dx = umin; 1991 else if (dx < 0.0 && dx > -umin) dx = -umin; 1992 #else 1993 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1994 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1995 #endif 1996 dx *= epsilon; 1997 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 1998 w3_array[col] += dx; 1999 } 2000 w3_array = w3_array + start; ierr = VecRestoreArrayFast(w3,&w3_array);CHKERRQ(ierr); 2001 2002 /* 2003 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2004 */ 2005 2006 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2007 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2008 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2009 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2010 2011 /* 2012 Loop over rows of vector, putting results into Jacobian matrix 2013 */ 2014 ierr = VecGetArrayFast(w2,&y);CHKERRQ(ierr); 2015 for (l=0; l<coloring->nrows[k]; l++) { 2016 row = coloring->rows[k][l]; 2017 col = coloring->columnsforrow[k][l]; 2018 y[row] *= vscale_array[vscaleforrow[k][l]]; 2019 srow = row + start; 2020 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2021 } 2022 ierr = VecRestoreArrayFast(w2,&y);CHKERRQ(ierr); 2023 } 2024 coloring->currentcolor = k; 2025 ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2026 xx = xx + start; ierr = VecRestoreArrayFast(x1,&xx);CHKERRQ(ierr); 2027 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2028 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2029 PetscFunctionReturn(0); 2030 } 2031 2032 #include "petscblaslapack.h" 2033 #undef __FUNCT__ 2034 #define __FUNCT__ "MatAXPY_SeqAIJ" 2035 int MatAXPY_SeqAIJ(const PetscScalar *a,Mat X,Mat Y,MatStructure str) 2036 { 2037 int ierr,one=1,i; 2038 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2039 2040 PetscFunctionBegin; 2041 if (str == SAME_NONZERO_PATTERN) { 2042 BLaxpy_(&x->nz,(PetscScalar*)a,x->a,&one,y->a,&one); 2043 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2044 if (y->xtoy && y->XtoY != X) { 2045 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2046 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2047 } 2048 if (!y->xtoy) { /* get xtoy */ 2049 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2050 y->XtoY = X; 2051 } 2052 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2053 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2054 } else { 2055 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2056 } 2057 PetscFunctionReturn(0); 2058 } 2059 2060 /* -------------------------------------------------------------------*/ 2061 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2062 MatGetRow_SeqAIJ, 2063 MatRestoreRow_SeqAIJ, 2064 MatMult_SeqAIJ, 2065 MatMultAdd_SeqAIJ, 2066 MatMultTranspose_SeqAIJ, 2067 MatMultTransposeAdd_SeqAIJ, 2068 MatSolve_SeqAIJ, 2069 MatSolveAdd_SeqAIJ, 2070 MatSolveTranspose_SeqAIJ, 2071 MatSolveTransposeAdd_SeqAIJ, 2072 MatLUFactor_SeqAIJ, 2073 0, 2074 MatRelax_SeqAIJ, 2075 MatTranspose_SeqAIJ, 2076 MatGetInfo_SeqAIJ, 2077 MatEqual_SeqAIJ, 2078 MatGetDiagonal_SeqAIJ, 2079 MatDiagonalScale_SeqAIJ, 2080 MatNorm_SeqAIJ, 2081 0, 2082 MatAssemblyEnd_SeqAIJ, 2083 MatCompress_SeqAIJ, 2084 MatSetOption_SeqAIJ, 2085 MatZeroEntries_SeqAIJ, 2086 MatZeroRows_SeqAIJ, 2087 MatLUFactorSymbolic_SeqAIJ, 2088 MatLUFactorNumeric_SeqAIJ, 2089 MatCholeskyFactorSymbolic_SeqAIJ, 2090 MatCholeskyFactorNumeric_SeqAIJ, 2091 MatSetUpPreallocation_SeqAIJ, 2092 MatILUFactorSymbolic_SeqAIJ, 2093 MatICCFactorSymbolic_SeqAIJ, 2094 MatGetArray_SeqAIJ, 2095 MatRestoreArray_SeqAIJ, 2096 MatDuplicate_SeqAIJ, 2097 0, 2098 0, 2099 MatILUFactor_SeqAIJ, 2100 0, 2101 MatAXPY_SeqAIJ, 2102 MatGetSubMatrices_SeqAIJ, 2103 MatIncreaseOverlap_SeqAIJ, 2104 MatGetValues_SeqAIJ, 2105 MatCopy_SeqAIJ, 2106 MatPrintHelp_SeqAIJ, 2107 MatScale_SeqAIJ, 2108 0, 2109 0, 2110 MatILUDTFactor_SeqAIJ, 2111 MatGetBlockSize_SeqAIJ, 2112 MatGetRowIJ_SeqAIJ, 2113 MatRestoreRowIJ_SeqAIJ, 2114 MatGetColumnIJ_SeqAIJ, 2115 MatRestoreColumnIJ_SeqAIJ, 2116 MatFDColoringCreate_SeqAIJ, 2117 0, 2118 0, 2119 MatPermute_SeqAIJ, 2120 0, 2121 0, 2122 MatDestroy_SeqAIJ, 2123 MatView_SeqAIJ, 2124 MatGetPetscMaps_Petsc, 2125 0, 2126 0, 2127 0, 2128 0, 2129 0, 2130 0, 2131 0, 2132 0, 2133 MatSetColoring_SeqAIJ, 2134 MatSetValuesAdic_SeqAIJ, 2135 MatSetValuesAdifor_SeqAIJ, 2136 MatFDColoringApply_SeqAIJ}; 2137 2138 EXTERN_C_BEGIN 2139 #undef __FUNCT__ 2140 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2141 2142 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2143 { 2144 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2145 int i,nz,n; 2146 2147 PetscFunctionBegin; 2148 2149 nz = aij->maxnz; 2150 n = mat->n; 2151 for (i=0; i<nz; i++) { 2152 aij->j[i] = indices[i]; 2153 } 2154 aij->nz = nz; 2155 for (i=0; i<n; i++) { 2156 aij->ilen[i] = aij->imax[i]; 2157 } 2158 2159 PetscFunctionReturn(0); 2160 } 2161 EXTERN_C_END 2162 2163 #undef __FUNCT__ 2164 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2165 /*@ 2166 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2167 in the matrix. 2168 2169 Input Parameters: 2170 + mat - the SeqAIJ matrix 2171 - indices - the column indices 2172 2173 Level: advanced 2174 2175 Notes: 2176 This can be called if you have precomputed the nonzero structure of the 2177 matrix and want to provide it to the matrix object to improve the performance 2178 of the MatSetValues() operation. 2179 2180 You MUST have set the correct numbers of nonzeros per row in the call to 2181 MatCreateSeqAIJ(). 2182 2183 MUST be called before any calls to MatSetValues(); 2184 2185 The indices should start with zero, not one. 2186 2187 @*/ 2188 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2189 { 2190 int ierr,(*f)(Mat,int *); 2191 2192 PetscFunctionBegin; 2193 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2194 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2195 if (f) { 2196 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2197 } else { 2198 SETERRQ(1,"Wrong type of matrix to set column indices"); 2199 } 2200 PetscFunctionReturn(0); 2201 } 2202 2203 /* ----------------------------------------------------------------------------------------*/ 2204 2205 EXTERN_C_BEGIN 2206 #undef __FUNCT__ 2207 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2208 int MatStoreValues_SeqAIJ(Mat mat) 2209 { 2210 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2211 size_t nz = aij->i[mat->m],ierr; 2212 2213 PetscFunctionBegin; 2214 if (aij->nonew != 1) { 2215 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2216 } 2217 2218 /* allocate space for values if not already there */ 2219 if (!aij->saved_values) { 2220 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2221 } 2222 2223 /* copy values over */ 2224 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2225 PetscFunctionReturn(0); 2226 } 2227 EXTERN_C_END 2228 2229 #undef __FUNCT__ 2230 #define __FUNCT__ "MatStoreValues" 2231 /*@ 2232 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2233 example, reuse of the linear part of a Jacobian, while recomputing the 2234 nonlinear portion. 2235 2236 Collect on Mat 2237 2238 Input Parameters: 2239 . mat - the matrix (currently on AIJ matrices support this option) 2240 2241 Level: advanced 2242 2243 Common Usage, with SNESSolve(): 2244 $ Create Jacobian matrix 2245 $ Set linear terms into matrix 2246 $ Apply boundary conditions to matrix, at this time matrix must have 2247 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2248 $ boundary conditions again will not change the nonzero structure 2249 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2250 $ ierr = MatStoreValues(mat); 2251 $ Call SNESSetJacobian() with matrix 2252 $ In your Jacobian routine 2253 $ ierr = MatRetrieveValues(mat); 2254 $ Set nonlinear terms in matrix 2255 2256 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2257 $ // build linear portion of Jacobian 2258 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2259 $ ierr = MatStoreValues(mat); 2260 $ loop over nonlinear iterations 2261 $ ierr = MatRetrieveValues(mat); 2262 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2263 $ // call MatAssemblyBegin/End() on matrix 2264 $ Solve linear system with Jacobian 2265 $ endloop 2266 2267 Notes: 2268 Matrix must already be assemblied before calling this routine 2269 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2270 calling this routine. 2271 2272 .seealso: MatRetrieveValues() 2273 2274 @*/ 2275 int MatStoreValues(Mat mat) 2276 { 2277 int ierr,(*f)(Mat); 2278 2279 PetscFunctionBegin; 2280 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2281 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2282 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2283 2284 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2285 if (f) { 2286 ierr = (*f)(mat);CHKERRQ(ierr); 2287 } else { 2288 SETERRQ(1,"Wrong type of matrix to store values"); 2289 } 2290 PetscFunctionReturn(0); 2291 } 2292 2293 EXTERN_C_BEGIN 2294 #undef __FUNCT__ 2295 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2296 int MatRetrieveValues_SeqAIJ(Mat mat) 2297 { 2298 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2299 int nz = aij->i[mat->m],ierr; 2300 2301 PetscFunctionBegin; 2302 if (aij->nonew != 1) { 2303 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2304 } 2305 if (!aij->saved_values) { 2306 SETERRQ(1,"Must call MatStoreValues(A);first"); 2307 } 2308 2309 /* copy values over */ 2310 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2311 PetscFunctionReturn(0); 2312 } 2313 EXTERN_C_END 2314 2315 #undef __FUNCT__ 2316 #define __FUNCT__ "MatRetrieveValues" 2317 /*@ 2318 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2319 example, reuse of the linear part of a Jacobian, while recomputing the 2320 nonlinear portion. 2321 2322 Collect on Mat 2323 2324 Input Parameters: 2325 . mat - the matrix (currently on AIJ matrices support this option) 2326 2327 Level: advanced 2328 2329 .seealso: MatStoreValues() 2330 2331 @*/ 2332 int MatRetrieveValues(Mat mat) 2333 { 2334 int ierr,(*f)(Mat); 2335 2336 PetscFunctionBegin; 2337 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2338 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2339 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2340 2341 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2342 if (f) { 2343 ierr = (*f)(mat);CHKERRQ(ierr); 2344 } else { 2345 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2346 } 2347 PetscFunctionReturn(0); 2348 } 2349 2350 /* 2351 This allows SeqAIJ matrices to be passed to the matlab engine 2352 */ 2353 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2354 #include "engine.h" /* Matlab include file */ 2355 #include "mex.h" /* Matlab include file */ 2356 EXTERN_C_BEGIN 2357 #undef __FUNCT__ 2358 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2359 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 2360 { 2361 int ierr; 2362 Mat B = (Mat)obj; 2363 mxArray *mat; 2364 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2365 2366 PetscFunctionBegin; 2367 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2368 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2369 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2370 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2371 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2372 2373 /* Matlab indices start at 0 for sparse (what a surprise) */ 2374 2375 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2376 mxSetName(mat,obj->name); 2377 engPutArray((Engine *)mengine,mat); 2378 PetscFunctionReturn(0); 2379 } 2380 EXTERN_C_END 2381 2382 EXTERN_C_BEGIN 2383 #undef __FUNCT__ 2384 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2385 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 2386 { 2387 int ierr,ii; 2388 Mat mat = (Mat)obj; 2389 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2390 mxArray *mmat; 2391 2392 PetscFunctionBegin; 2393 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2394 2395 mmat = engGetArray((Engine *)mengine,obj->name); 2396 2397 aij->nz = (mxGetJc(mmat))[mat->m]; 2398 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2399 aij->j = (int*)(aij->a + aij->nz); 2400 aij->i = aij->j + aij->nz; 2401 aij->singlemalloc = PETSC_TRUE; 2402 aij->freedata = PETSC_TRUE; 2403 2404 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2405 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2406 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2407 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2408 2409 for (ii=0; ii<mat->m; ii++) { 2410 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2411 } 2412 2413 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2414 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2415 2416 PetscFunctionReturn(0); 2417 } 2418 EXTERN_C_END 2419 #endif 2420 2421 /* --------------------------------------------------------------------------------*/ 2422 #undef __FUNCT__ 2423 #define __FUNCT__ "MatCreateSeqAIJ" 2424 /*@C 2425 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2426 (the default parallel PETSc format). For good matrix assembly performance 2427 the user should preallocate the matrix storage by setting the parameter nz 2428 (or the array nnz). By setting these parameters accurately, performance 2429 during matrix assembly can be increased by more than a factor of 50. 2430 2431 Collective on MPI_Comm 2432 2433 Input Parameters: 2434 + comm - MPI communicator, set to PETSC_COMM_SELF 2435 . m - number of rows 2436 . n - number of columns 2437 . nz - number of nonzeros per row (same for all rows) 2438 - nnz - array containing the number of nonzeros in the various rows 2439 (possibly different for each row) or PETSC_NULL 2440 2441 Output Parameter: 2442 . A - the matrix 2443 2444 Notes: 2445 The AIJ format (also called the Yale sparse matrix format or 2446 compressed row storage), is fully compatible with standard Fortran 77 2447 storage. That is, the stored row and column indices can begin at 2448 either one (as in Fortran) or zero. See the users' manual for details. 2449 2450 Specify the preallocated storage with either nz or nnz (not both). 2451 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2452 allocation. For large problems you MUST preallocate memory or you 2453 will get TERRIBLE performance, see the users' manual chapter on matrices. 2454 2455 By default, this format uses inodes (identical nodes) when possible, to 2456 improve numerical efficiency of matrix-vector products and solves. We 2457 search for consecutive rows with the same nonzero structure, thereby 2458 reusing matrix information to achieve increased efficiency. 2459 2460 Options Database Keys: 2461 + -mat_aij_no_inode - Do not use inodes 2462 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2463 - -mat_aij_oneindex - Internally use indexing starting at 1 2464 rather than 0. Note that when calling MatSetValues(), 2465 the user still MUST index entries starting at 0! 2466 2467 Level: intermediate 2468 2469 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2470 2471 @*/ 2472 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,const int nnz[],Mat *A) 2473 { 2474 int ierr; 2475 2476 PetscFunctionBegin; 2477 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2478 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2479 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2480 PetscFunctionReturn(0); 2481 } 2482 2483 #define SKIP_ALLOCATION -4 2484 2485 #undef __FUNCT__ 2486 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2487 /*@C 2488 MatSeqAIJSetPreallocation - For good matrix assembly performance 2489 the user should preallocate the matrix storage by setting the parameter nz 2490 (or the array nnz). By setting these parameters accurately, performance 2491 during matrix assembly can be increased by more than a factor of 50. 2492 2493 Collective on MPI_Comm 2494 2495 Input Parameters: 2496 + comm - MPI communicator, set to PETSC_COMM_SELF 2497 . m - number of rows 2498 . n - number of columns 2499 . nz - number of nonzeros per row (same for all rows) 2500 - nnz - array containing the number of nonzeros in the various rows 2501 (possibly different for each row) or PETSC_NULL 2502 2503 Output Parameter: 2504 . A - the matrix 2505 2506 Notes: 2507 The AIJ format (also called the Yale sparse matrix format or 2508 compressed row storage), is fully compatible with standard Fortran 77 2509 storage. That is, the stored row and column indices can begin at 2510 either one (as in Fortran) or zero. See the users' manual for details. 2511 2512 Specify the preallocated storage with either nz or nnz (not both). 2513 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2514 allocation. For large problems you MUST preallocate memory or you 2515 will get TERRIBLE performance, see the users' manual chapter on matrices. 2516 2517 By default, this format uses inodes (identical nodes) when possible, to 2518 improve numerical efficiency of matrix-vector products and solves. We 2519 search for consecutive rows with the same nonzero structure, thereby 2520 reusing matrix information to achieve increased efficiency. 2521 2522 Options Database Keys: 2523 + -mat_aij_no_inode - Do not use inodes 2524 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2525 - -mat_aij_oneindex - Internally use indexing starting at 1 2526 rather than 0. Note that when calling MatSetValues(), 2527 the user still MUST index entries starting at 0! 2528 2529 Level: intermediate 2530 2531 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2532 2533 @*/ 2534 int MatSeqAIJSetPreallocation(Mat B,int nz,const int nnz[]) 2535 { 2536 int ierr,(*f)(Mat,int,const int[]); 2537 2538 PetscFunctionBegin; 2539 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2540 if (f) { 2541 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2542 } 2543 PetscFunctionReturn(0); 2544 } 2545 2546 EXTERN_C_BEGIN 2547 #undef __FUNCT__ 2548 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2549 int MatSeqAIJSetPreallocation_SeqAIJ(Mat B,int nz,int *nnz) 2550 { 2551 Mat_SeqAIJ *b; 2552 size_t len = 0; 2553 PetscTruth skipallocation = PETSC_FALSE; 2554 int i,ierr; 2555 2556 PetscFunctionBegin; 2557 2558 if (nz == SKIP_ALLOCATION) { 2559 skipallocation = PETSC_TRUE; 2560 nz = 0; 2561 } 2562 2563 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2564 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2565 if (nnz) { 2566 for (i=0; i<B->m; i++) { 2567 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2568 if (nnz[i] > B->n) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->n); 2569 } 2570 } 2571 2572 B->preallocated = PETSC_TRUE; 2573 b = (Mat_SeqAIJ*)B->data; 2574 2575 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2576 if (!nnz) { 2577 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2578 else if (nz <= 0) nz = 1; 2579 for (i=0; i<B->m; i++) b->imax[i] = nz; 2580 nz = nz*B->m; 2581 } else { 2582 nz = 0; 2583 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2584 } 2585 2586 if (!skipallocation) { 2587 /* allocate the matrix space */ 2588 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2589 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2590 b->j = (int*)(b->a + nz); 2591 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2592 b->i = b->j + nz; 2593 b->i[0] = 0; 2594 for (i=1; i<B->m+1; i++) { 2595 b->i[i] = b->i[i-1] + b->imax[i-1]; 2596 } 2597 b->singlemalloc = PETSC_TRUE; 2598 b->freedata = PETSC_TRUE; 2599 } else { 2600 b->freedata = PETSC_FALSE; 2601 } 2602 2603 /* b->ilen will count nonzeros in each row so far. */ 2604 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2605 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2606 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2607 2608 b->nz = 0; 2609 b->maxnz = nz; 2610 B->info.nz_unneeded = (double)b->maxnz; 2611 PetscFunctionReturn(0); 2612 } 2613 EXTERN_C_END 2614 2615 EXTERN int RegisterApplyPtAPRoutines_Private(Mat); 2616 2617 EXTERN_C_BEGIN 2618 EXTERN int MatConvert_SeqAIJ_SeqSBAIJ(Mat,MatType,Mat*); 2619 EXTERN int MatConvert_SeqAIJ_SeqBAIJ(Mat,MatType,Mat*); 2620 EXTERN int MatReorderForNonzeroDiagonal_SeqAIJ(Mat,PetscReal,IS,IS); 2621 EXTERN int MatAdjustForInodes_SeqAIJ(Mat,IS*,IS*); 2622 EXTERN int MatSeqAIJGetInodeSizes_SeqAIJ(Mat,int*,int*[],int*); 2623 EXTERN_C_END 2624 2625 EXTERN_C_BEGIN 2626 #undef __FUNCT__ 2627 #define __FUNCT__ "MatCreate_SeqAIJ" 2628 int MatCreate_SeqAIJ(Mat B) 2629 { 2630 Mat_SeqAIJ *b; 2631 int ierr,size; 2632 PetscTruth flg; 2633 2634 PetscFunctionBegin; 2635 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2636 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2637 2638 B->m = B->M = PetscMax(B->m,B->M); 2639 B->n = B->N = PetscMax(B->n,B->N); 2640 2641 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2642 B->data = (void*)b; 2643 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2644 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2645 B->factor = 0; 2646 B->lupivotthreshold = 1.0; 2647 B->mapping = 0; 2648 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2649 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2650 b->row = 0; 2651 b->col = 0; 2652 b->icol = 0; 2653 b->reallocs = 0; 2654 2655 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2656 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2657 2658 b->sorted = PETSC_FALSE; 2659 b->ignorezeroentries = PETSC_FALSE; 2660 b->roworiented = PETSC_TRUE; 2661 b->nonew = 0; 2662 b->diag = 0; 2663 b->solve_work = 0; 2664 B->spptr = 0; 2665 b->inode.use = PETSC_TRUE; 2666 b->inode.node_count = 0; 2667 b->inode.size = 0; 2668 b->inode.limit = 5; 2669 b->inode.max_limit = 5; 2670 b->saved_values = 0; 2671 b->idiag = 0; 2672 b->ssor = 0; 2673 b->keepzeroedrows = PETSC_FALSE; 2674 b->xtoy = 0; 2675 b->XtoY = 0; 2676 2677 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2678 2679 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2680 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2681 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2682 "MatSeqAIJSetColumnIndices_SeqAIJ", 2683 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2684 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2685 "MatStoreValues_SeqAIJ", 2686 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2687 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2688 "MatRetrieveValues_SeqAIJ", 2689 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2690 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2691 "MatConvert_SeqAIJ_SeqSBAIJ", 2692 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2693 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2694 "MatConvert_SeqAIJ_SeqBAIJ", 2695 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2696 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsSymmetric_C", 2697 "MatIsSymmetric_SeqAIJ", 2698 MatIsSymmetric_SeqAIJ);CHKERRQ(ierr); 2699 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2700 "MatSeqAIJSetPreallocation_SeqAIJ", 2701 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2702 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2703 "MatReorderForNonzeroDiagonal_SeqAIJ", 2704 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2705 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2706 "MatAdjustForInodes_SeqAIJ", 2707 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2708 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2709 "MatSeqAIJGetInodeSizes_SeqAIJ", 2710 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2711 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2712 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2713 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2714 #endif 2715 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2716 PetscFunctionReturn(0); 2717 } 2718 EXTERN_C_END 2719 2720 #undef __FUNCT__ 2721 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2722 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2723 { 2724 Mat C; 2725 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2726 int i,m = A->m,ierr; 2727 size_t len; 2728 2729 PetscFunctionBegin; 2730 *B = 0; 2731 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2732 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2733 c = (Mat_SeqAIJ*)C->data; 2734 2735 C->factor = A->factor; 2736 c->row = 0; 2737 c->col = 0; 2738 c->icol = 0; 2739 c->keepzeroedrows = a->keepzeroedrows; 2740 C->assembled = PETSC_TRUE; 2741 2742 C->M = A->m; 2743 C->N = A->n; 2744 2745 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2746 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2747 for (i=0; i<m; i++) { 2748 c->imax[i] = a->imax[i]; 2749 c->ilen[i] = a->ilen[i]; 2750 } 2751 2752 /* allocate the matrix space */ 2753 c->singlemalloc = PETSC_TRUE; 2754 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2755 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2756 c->j = (int*)(c->a + a->i[m] ); 2757 c->i = c->j + a->i[m]; 2758 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2759 if (m > 0) { 2760 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(int));CHKERRQ(ierr); 2761 if (cpvalues == MAT_COPY_VALUES) { 2762 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2763 } else { 2764 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2765 } 2766 } 2767 2768 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2769 c->sorted = a->sorted; 2770 c->roworiented = a->roworiented; 2771 c->nonew = a->nonew; 2772 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2773 c->saved_values = 0; 2774 c->idiag = 0; 2775 c->ssor = 0; 2776 c->ignorezeroentries = a->ignorezeroentries; 2777 c->freedata = PETSC_TRUE; 2778 2779 if (a->diag) { 2780 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2781 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2782 for (i=0; i<m; i++) { 2783 c->diag[i] = a->diag[i]; 2784 } 2785 } else c->diag = 0; 2786 c->inode.use = a->inode.use; 2787 c->inode.limit = a->inode.limit; 2788 c->inode.max_limit = a->inode.max_limit; 2789 if (a->inode.size){ 2790 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2791 c->inode.node_count = a->inode.node_count; 2792 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2793 } else { 2794 c->inode.size = 0; 2795 c->inode.node_count = 0; 2796 } 2797 c->nz = a->nz; 2798 c->maxnz = a->maxnz; 2799 c->solve_work = 0; 2800 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2801 C->preallocated = PETSC_TRUE; 2802 2803 *B = C; 2804 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2805 PetscFunctionReturn(0); 2806 } 2807 2808 EXTERN_C_BEGIN 2809 #undef __FUNCT__ 2810 #define __FUNCT__ "MatLoad_SeqAIJ" 2811 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2812 { 2813 Mat_SeqAIJ *a; 2814 Mat B; 2815 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N; 2816 MPI_Comm comm; 2817 2818 PetscFunctionBegin; 2819 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2820 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2821 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2822 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2823 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2824 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2825 M = header[1]; N = header[2]; nz = header[3]; 2826 2827 if (nz < 0) { 2828 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2829 } 2830 2831 /* read in row lengths */ 2832 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2833 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2834 2835 /* create our matrix */ 2836 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2837 ierr = MatSetType(B,type);CHKERRQ(ierr); 2838 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2839 a = (Mat_SeqAIJ*)B->data; 2840 2841 /* read in column indices and adjust for Fortran indexing*/ 2842 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2843 2844 /* read in nonzero values */ 2845 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2846 2847 /* set matrix "i" values */ 2848 a->i[0] = 0; 2849 for (i=1; i<= M; i++) { 2850 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2851 a->ilen[i-1] = rowlengths[i-1]; 2852 } 2853 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2854 2855 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2856 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2857 *A = B; 2858 PetscFunctionReturn(0); 2859 } 2860 EXTERN_C_END 2861 2862 #undef __FUNCT__ 2863 #define __FUNCT__ "MatEqual_SeqAIJ" 2864 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2865 { 2866 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2867 int ierr; 2868 2869 PetscFunctionBegin; 2870 /* If the matrix dimensions are not equal,or no of nonzeros */ 2871 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2872 *flg = PETSC_FALSE; 2873 PetscFunctionReturn(0); 2874 } 2875 2876 /* if the a->i are the same */ 2877 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2878 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2879 2880 /* if a->j are the same */ 2881 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2882 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2883 2884 /* if a->a are the same */ 2885 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2886 2887 PetscFunctionReturn(0); 2888 2889 } 2890 2891 #undef __FUNCT__ 2892 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2893 /*@C 2894 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2895 provided by the user. 2896 2897 Coolective on MPI_Comm 2898 2899 Input Parameters: 2900 + comm - must be an MPI communicator of size 1 2901 . m - number of rows 2902 . n - number of columns 2903 . i - row indices 2904 . j - column indices 2905 - a - matrix values 2906 2907 Output Parameter: 2908 . mat - the matrix 2909 2910 Level: intermediate 2911 2912 Notes: 2913 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2914 once the matrix is destroyed 2915 2916 You cannot set new nonzero locations into this matrix, that will generate an error. 2917 2918 The i and j indices are 0 based 2919 2920 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2921 2922 @*/ 2923 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2924 { 2925 int ierr,ii; 2926 Mat_SeqAIJ *aij; 2927 2928 PetscFunctionBegin; 2929 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2930 aij = (Mat_SeqAIJ*)(*mat)->data; 2931 2932 if (i[0] != 0) { 2933 SETERRQ(1,"i (row indices) must start with 0"); 2934 } 2935 aij->i = i; 2936 aij->j = j; 2937 aij->a = a; 2938 aij->singlemalloc = PETSC_FALSE; 2939 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2940 aij->freedata = PETSC_FALSE; 2941 2942 for (ii=0; ii<m; ii++) { 2943 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2944 #if defined(PETSC_USE_BOPT_g) 2945 if (i[ii+1] - i[ii] < 0) SETERRQ2(1,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2946 #endif 2947 } 2948 #if defined(PETSC_USE_BOPT_g) 2949 for (ii=0; ii<aij->i[m]; ii++) { 2950 if (j[ii] < 0) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2951 if (j[ii] > n - 1) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2952 } 2953 #endif 2954 2955 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2956 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2957 PetscFunctionReturn(0); 2958 } 2959 2960 #undef __FUNCT__ 2961 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2962 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2963 { 2964 int ierr; 2965 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2966 2967 PetscFunctionBegin; 2968 if (coloring->ctype == IS_COLORING_LOCAL) { 2969 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2970 a->coloring = coloring; 2971 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 2972 int i,*larray; 2973 ISColoring ocoloring; 2974 ISColoringValue *colors; 2975 2976 /* set coloring for diagonal portion */ 2977 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 2978 for (i=0; i<A->n; i++) { 2979 larray[i] = i; 2980 } 2981 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 2982 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 2983 for (i=0; i<A->n; i++) { 2984 colors[i] = coloring->colors[larray[i]]; 2985 } 2986 ierr = PetscFree(larray);CHKERRQ(ierr); 2987 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 2988 a->coloring = ocoloring; 2989 } 2990 PetscFunctionReturn(0); 2991 } 2992 2993 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2994 EXTERN_C_BEGIN 2995 #include "adic/ad_utils.h" 2996 EXTERN_C_END 2997 2998 #undef __FUNCT__ 2999 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3000 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3001 { 3002 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3003 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3004 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3005 ISColoringValue *color; 3006 3007 PetscFunctionBegin; 3008 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3009 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3010 color = a->coloring->colors; 3011 /* loop over rows */ 3012 for (i=0; i<m; i++) { 3013 nz = ii[i+1] - ii[i]; 3014 /* loop over columns putting computed value into matrix */ 3015 for (j=0; j<nz; j++) { 3016 *v++ = values[color[*jj++]]; 3017 } 3018 values += nlen; /* jump to next row of derivatives */ 3019 } 3020 PetscFunctionReturn(0); 3021 } 3022 3023 #else 3024 3025 #undef __FUNCT__ 3026 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3027 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3028 { 3029 PetscFunctionBegin; 3030 SETERRQ(1,"PETSc installed without ADIC"); 3031 } 3032 3033 #endif 3034 3035 #undef __FUNCT__ 3036 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3037 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3038 { 3039 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3040 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3041 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3042 ISColoringValue *color; 3043 3044 PetscFunctionBegin; 3045 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3046 color = a->coloring->colors; 3047 /* loop over rows */ 3048 for (i=0; i<m; i++) { 3049 nz = ii[i+1] - ii[i]; 3050 /* loop over columns putting computed value into matrix */ 3051 for (j=0; j<nz; j++) { 3052 *v++ = values[color[*jj++]]; 3053 } 3054 values += nl; /* jump to next row of derivatives */ 3055 } 3056 PetscFunctionReturn(0); 3057 } 3058 3059