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