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