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,*b,*bs, d,*xs,sum,*v = a->a,*t=0,scale,*ts,*xb,*idiag=0,*mdiag; 1027 int ierr,*idx,*diag,n = A->n,m = A->m,i; 1028 1029 PetscFunctionBegin; 1030 its = its*lits; 1031 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits); 1032 1033 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1034 diag = a->diag; 1035 if (!a->idiag) { 1036 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1037 a->ssor = a->idiag + m; 1038 mdiag = a->ssor + m; 1039 1040 v = a->a; 1041 1042 /* this is wrong when fshift omega changes each iteration */ 1043 if (omega == 1.0 && fshift == 0.0) { 1044 for (i=0; i<m; i++) { 1045 mdiag[i] = v[diag[i]]; 1046 a->idiag[i] = 1.0/v[diag[i]]; 1047 } 1048 } else { 1049 for (i=0; i<m; i++) { 1050 mdiag[i] = v[diag[i]]; 1051 a->idiag[i] = omega/(fshift + v[diag[i]]);} 1052 } 1053 1054 } 1055 t = a->ssor; 1056 idiag = a->idiag; 1057 mdiag = a->idiag + 2*m; 1058 1059 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 1060 if (xx != bb) { 1061 ierr = VecGetArrayFast(bb,&b);CHKERRQ(ierr); 1062 } else { 1063 b = x; 1064 } 1065 1066 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1067 xs = x; 1068 if (flag == SOR_APPLY_UPPER) { 1069 /* apply (U + D/omega) to the vector */ 1070 bs = b; 1071 for (i=0; i<m; i++) { 1072 d = fshift + a->a[diag[i]]; 1073 n = a->i[i+1] - diag[i] - 1; 1074 idx = a->j + diag[i] + 1; 1075 v = a->a + diag[i] + 1; 1076 sum = b[i]*d/omega; 1077 SPARSEDENSEDOT(sum,bs,v,idx,n); 1078 x[i] = sum; 1079 } 1080 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1081 if (bb != xx) {ierr = VecRestoreArrayFast(bb,&b);CHKERRQ(ierr);} 1082 PetscLogFlops(a->nz); 1083 PetscFunctionReturn(0); 1084 } 1085 1086 1087 /* Let A = L + U + D; where L is lower trianglar, 1088 U is upper triangular, E is diagonal; This routine applies 1089 1090 (L + E)^{-1} A (U + E)^{-1} 1091 1092 to a vector efficiently using Eisenstat's trick. This is for 1093 the case of SSOR preconditioner, so E is D/omega where omega 1094 is the relaxation factor. 1095 */ 1096 1097 if (flag == SOR_APPLY_LOWER) { 1098 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1099 } else if (flag & SOR_EISENSTAT) { 1100 /* Let A = L + U + D; where L is lower trianglar, 1101 U is upper triangular, E is diagonal; This routine applies 1102 1103 (L + E)^{-1} A (U + E)^{-1} 1104 1105 to a vector efficiently using Eisenstat's trick. This is for 1106 the case of SSOR preconditioner, so E is D/omega where omega 1107 is the relaxation factor. 1108 */ 1109 scale = (2.0/omega) - 1.0; 1110 1111 /* x = (E + U)^{-1} b */ 1112 for (i=m-1; i>=0; i--) { 1113 n = a->i[i+1] - diag[i] - 1; 1114 idx = a->j + diag[i] + 1; 1115 v = a->a + diag[i] + 1; 1116 sum = b[i]; 1117 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1118 x[i] = sum*idiag[i]; 1119 } 1120 1121 /* t = b - (2*E - D)x */ 1122 v = a->a; 1123 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1124 1125 /* t = (E + L)^{-1}t */ 1126 ts = t; 1127 diag = a->diag; 1128 for (i=0; i<m; i++) { 1129 n = diag[i] - a->i[i]; 1130 idx = a->j + a->i[i]; 1131 v = a->a + a->i[i]; 1132 sum = t[i]; 1133 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1134 t[i] = sum*idiag[i]; 1135 /* x = x + t */ 1136 x[i] += t[i]; 1137 } 1138 1139 PetscLogFlops(6*m-1 + 2*a->nz); 1140 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1141 if (bb != xx) {ierr = VecRestoreArrayFast(bb,&b);CHKERRQ(ierr);} 1142 PetscFunctionReturn(0); 1143 } 1144 if (flag & SOR_ZERO_INITIAL_GUESS) { 1145 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1146 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1147 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1148 #else 1149 for (i=0; i<m; i++) { 1150 n = diag[i] - a->i[i]; 1151 idx = a->j + a->i[i]; 1152 v = a->a + a->i[i]; 1153 sum = b[i]; 1154 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1155 x[i] = sum*idiag[i]; 1156 } 1157 #endif 1158 xb = x; 1159 PetscLogFlops(a->nz); 1160 } else xb = b; 1161 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1162 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1163 for (i=0; i<m; i++) { 1164 x[i] *= mdiag[i]; 1165 } 1166 PetscLogFlops(m); 1167 } 1168 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1169 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1170 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,diag,a->a,xb); 1171 #else 1172 for (i=m-1; i>=0; i--) { 1173 n = a->i[i+1] - diag[i] - 1; 1174 idx = a->j + diag[i] + 1; 1175 v = a->a + diag[i] + 1; 1176 sum = xb[i]; 1177 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1178 x[i] = sum*idiag[i]; 1179 } 1180 #endif 1181 PetscLogFlops(a->nz); 1182 } 1183 its--; 1184 } 1185 while (its--) { 1186 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1187 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1188 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1189 #else 1190 for (i=0; i<m; i++) { 1191 d = fshift + a->a[diag[i]]; 1192 n = a->i[i+1] - a->i[i]; 1193 idx = a->j + a->i[i]; 1194 v = a->a + a->i[i]; 1195 sum = b[i]; 1196 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1197 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1198 } 1199 #endif 1200 PetscLogFlops(a->nz); 1201 } 1202 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1203 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1204 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1205 #else 1206 for (i=m-1; i>=0; i--) { 1207 d = fshift + a->a[diag[i]]; 1208 n = a->i[i+1] - a->i[i]; 1209 idx = a->j + a->i[i]; 1210 v = a->a + a->i[i]; 1211 sum = b[i]; 1212 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1213 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1214 } 1215 #endif 1216 PetscLogFlops(a->nz); 1217 } 1218 } 1219 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1220 if (bb != xx) {ierr = VecRestoreArrayFast(bb,&b);CHKERRQ(ierr);} 1221 PetscFunctionReturn(0); 1222 } 1223 1224 #undef __FUNCT__ 1225 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1226 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1227 { 1228 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1229 1230 PetscFunctionBegin; 1231 info->rows_global = (double)A->m; 1232 info->columns_global = (double)A->n; 1233 info->rows_local = (double)A->m; 1234 info->columns_local = (double)A->n; 1235 info->block_size = 1.0; 1236 info->nz_allocated = (double)a->maxnz; 1237 info->nz_used = (double)a->nz; 1238 info->nz_unneeded = (double)(a->maxnz - a->nz); 1239 info->assemblies = (double)A->num_ass; 1240 info->mallocs = (double)a->reallocs; 1241 info->memory = A->mem; 1242 if (A->factor) { 1243 info->fill_ratio_given = A->info.fill_ratio_given; 1244 info->fill_ratio_needed = A->info.fill_ratio_needed; 1245 info->factor_mallocs = A->info.factor_mallocs; 1246 } else { 1247 info->fill_ratio_given = 0; 1248 info->fill_ratio_needed = 0; 1249 info->factor_mallocs = 0; 1250 } 1251 PetscFunctionReturn(0); 1252 } 1253 1254 EXTERN int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,MatFactorInfo*,Mat*); 1255 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1256 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatFactorInfo*); 1257 EXTERN int MatSolve_SeqAIJ(Mat,Vec,Vec); 1258 EXTERN int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1259 EXTERN int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1260 EXTERN int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1261 1262 #undef __FUNCT__ 1263 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1264 int MatZeroRows_SeqAIJ(Mat A,IS is,const PetscScalar *diag) 1265 { 1266 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1267 int i,ierr,N,*rows,m = A->m - 1; 1268 1269 PetscFunctionBegin; 1270 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1271 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1272 if (a->keepzeroedrows) { 1273 for (i=0; i<N; i++) { 1274 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1275 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1276 } 1277 if (diag) { 1278 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1279 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1280 for (i=0; i<N; i++) { 1281 a->a[a->diag[rows[i]]] = *diag; 1282 } 1283 } 1284 } else { 1285 if (diag) { 1286 for (i=0; i<N; i++) { 1287 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1288 if (a->ilen[rows[i]] > 0) { 1289 a->ilen[rows[i]] = 1; 1290 a->a[a->i[rows[i]]] = *diag; 1291 a->j[a->i[rows[i]]] = rows[i]; 1292 } else { /* in case row was completely empty */ 1293 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1294 } 1295 } 1296 } else { 1297 for (i=0; i<N; i++) { 1298 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1299 a->ilen[rows[i]] = 0; 1300 } 1301 } 1302 } 1303 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1304 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1305 PetscFunctionReturn(0); 1306 } 1307 1308 #undef __FUNCT__ 1309 #define __FUNCT__ "MatGetRow_SeqAIJ" 1310 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1311 { 1312 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1313 int *itmp; 1314 1315 PetscFunctionBegin; 1316 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1317 1318 *nz = a->i[row+1] - a->i[row]; 1319 if (v) *v = a->a + a->i[row]; 1320 if (idx) { 1321 itmp = a->j + a->i[row]; 1322 if (*nz) { 1323 *idx = itmp; 1324 } 1325 else *idx = 0; 1326 } 1327 PetscFunctionReturn(0); 1328 } 1329 1330 /* remove this function? */ 1331 #undef __FUNCT__ 1332 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1333 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1334 { 1335 /* Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1336 int ierr; */ 1337 1338 PetscFunctionBegin; 1339 /* if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} */ 1340 PetscFunctionReturn(0); 1341 } 1342 1343 #undef __FUNCT__ 1344 #define __FUNCT__ "MatNorm_SeqAIJ" 1345 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1346 { 1347 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1348 PetscScalar *v = a->a; 1349 PetscReal sum = 0.0; 1350 int i,j,ierr; 1351 1352 PetscFunctionBegin; 1353 if (type == NORM_FROBENIUS) { 1354 for (i=0; i<a->nz; i++) { 1355 #if defined(PETSC_USE_COMPLEX) 1356 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1357 #else 1358 sum += (*v)*(*v); v++; 1359 #endif 1360 } 1361 *nrm = sqrt(sum); 1362 } else if (type == NORM_1) { 1363 PetscReal *tmp; 1364 int *jj = a->j; 1365 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1366 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1367 *nrm = 0.0; 1368 for (j=0; j<a->nz; j++) { 1369 tmp[*jj++] += PetscAbsScalar(*v); v++; 1370 } 1371 for (j=0; j<A->n; j++) { 1372 if (tmp[j] > *nrm) *nrm = tmp[j]; 1373 } 1374 ierr = PetscFree(tmp);CHKERRQ(ierr); 1375 } else if (type == NORM_INFINITY) { 1376 *nrm = 0.0; 1377 for (j=0; j<A->m; j++) { 1378 v = a->a + a->i[j]; 1379 sum = 0.0; 1380 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1381 sum += PetscAbsScalar(*v); v++; 1382 } 1383 if (sum > *nrm) *nrm = sum; 1384 } 1385 } else { 1386 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1387 } 1388 PetscFunctionReturn(0); 1389 } 1390 1391 #undef __FUNCT__ 1392 #define __FUNCT__ "MatTranspose_SeqAIJ" 1393 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1394 { 1395 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1396 Mat C; 1397 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1398 PetscScalar *array = a->a; 1399 1400 PetscFunctionBegin; 1401 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1402 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1403 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1404 1405 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1406 ierr = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1407 ierr = PetscFree(col);CHKERRQ(ierr); 1408 for (i=0; i<m; i++) { 1409 len = ai[i+1]-ai[i]; 1410 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1411 array += len; 1412 aj += len; 1413 } 1414 1415 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1416 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1417 1418 if (B) { 1419 *B = C; 1420 } else { 1421 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1422 } 1423 PetscFunctionReturn(0); 1424 } 1425 1426 EXTERN_C_BEGIN 1427 #undef __FUNCT__ 1428 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1429 int MatIsSymmetric_SeqAIJ(Mat A,Mat B,PetscTruth *f) 1430 { 1431 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1432 int *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1433 int ma,na,mb,nb, i,ierr; 1434 1435 PetscFunctionBegin; 1436 bij = (Mat_SeqAIJ *) B->data; 1437 1438 ierr = MatGetSize(A,&ma,&na); CHKERRQ(ierr); 1439 ierr = MatGetSize(B,&mb,&nb); CHKERRQ(ierr); 1440 if (ma!=nb || na!=mb) 1441 SETERRQ(1,"Incompatible A/B sizes for symmetry test"); 1442 aii = aij->i; bii = bij->i; 1443 adx = aij->j; bdx = bij->j; 1444 va = aij->a; vb = bij->a; 1445 ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr); 1446 ierr = PetscMalloc(mb*sizeof(int),&bptr); CHKERRQ(ierr); 1447 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1448 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1449 1450 *f = PETSC_TRUE; 1451 for (i=0; i<ma; i++) { 1452 /*printf("row %d spans %d--%d; we start @ %d\n", 1453 i,idx[ii[i]],idx[ii[i+1]-1],idx[aptr[i]]);*/ 1454 while (aptr[i]<aii[i+1]) { 1455 int idc,idr; PetscScalar vc,vr; 1456 /* column/row index/value */ 1457 idc = adx[aptr[i]]; idr = bdx[bptr[idc]]; 1458 vc = va[aptr[i]]; vr = vb[bptr[idc]]; 1459 /*printf("comparing %d: (%d,%d)=%e to (%d,%d)=%e\n", 1460 aptr[i],i,idc,vc,idc,idr,vr);*/ 1461 if (i!=idr || vc!=vr) { 1462 *f = PETSC_FALSE; goto done; 1463 } else { 1464 aptr[i]++; if (B || i!=idc) bptr[idc]++; 1465 } 1466 } 1467 } 1468 done: 1469 ierr = PetscFree(aptr); CHKERRQ(ierr); 1470 if (B) { 1471 ierr = PetscFree(bptr); CHKERRQ(ierr); 1472 } 1473 1474 PetscFunctionReturn(0); 1475 } 1476 EXTERN_C_END 1477 1478 #undef __FUNCT__ 1479 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1480 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1481 { 1482 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1483 PetscScalar *l,*r,x,*v; 1484 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1485 1486 PetscFunctionBegin; 1487 if (ll) { 1488 /* The local size is used so that VecMPI can be passed to this routine 1489 by MatDiagonalScale_MPIAIJ */ 1490 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1491 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1492 ierr = VecGetArrayFast(ll,&l);CHKERRQ(ierr); 1493 v = a->a; 1494 for (i=0; i<m; i++) { 1495 x = l[i]; 1496 M = a->i[i+1] - a->i[i]; 1497 for (j=0; j<M; j++) { (*v++) *= x;} 1498 } 1499 ierr = VecRestoreArrayFast(ll,&l);CHKERRQ(ierr); 1500 PetscLogFlops(nz); 1501 } 1502 if (rr) { 1503 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1504 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1505 ierr = VecGetArrayFast(rr,&r);CHKERRQ(ierr); 1506 v = a->a; jj = a->j; 1507 for (i=0; i<nz; i++) { 1508 (*v++) *= r[*jj++]; 1509 } 1510 ierr = VecRestoreArrayFast(rr,&r);CHKERRQ(ierr); 1511 PetscLogFlops(nz); 1512 } 1513 PetscFunctionReturn(0); 1514 } 1515 1516 #undef __FUNCT__ 1517 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1518 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1519 { 1520 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1521 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1522 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1523 int *irow,*icol,nrows,ncols; 1524 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1525 PetscScalar *a_new,*mat_a; 1526 Mat C; 1527 PetscTruth stride; 1528 1529 PetscFunctionBegin; 1530 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1531 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1532 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1533 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1534 1535 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1536 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1537 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1538 1539 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1540 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1541 if (stride && step == 1) { 1542 /* special case of contiguous rows */ 1543 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1544 starts = lens + nrows; 1545 /* loop over new rows determining lens and starting points */ 1546 for (i=0; i<nrows; i++) { 1547 kstart = ai[irow[i]]; 1548 kend = kstart + ailen[irow[i]]; 1549 for (k=kstart; k<kend; k++) { 1550 if (aj[k] >= first) { 1551 starts[i] = k; 1552 break; 1553 } 1554 } 1555 sum = 0; 1556 while (k < kend) { 1557 if (aj[k++] >= first+ncols) break; 1558 sum++; 1559 } 1560 lens[i] = sum; 1561 } 1562 /* create submatrix */ 1563 if (scall == MAT_REUSE_MATRIX) { 1564 int n_cols,n_rows; 1565 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1566 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1567 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1568 C = *B; 1569 } else { 1570 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1571 } 1572 c = (Mat_SeqAIJ*)C->data; 1573 1574 /* loop over rows inserting into submatrix */ 1575 a_new = c->a; 1576 j_new = c->j; 1577 i_new = c->i; 1578 1579 for (i=0; i<nrows; i++) { 1580 ii = starts[i]; 1581 lensi = lens[i]; 1582 for (k=0; k<lensi; k++) { 1583 *j_new++ = aj[ii+k] - first; 1584 } 1585 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1586 a_new += lensi; 1587 i_new[i+1] = i_new[i] + lensi; 1588 c->ilen[i] = lensi; 1589 } 1590 ierr = PetscFree(lens);CHKERRQ(ierr); 1591 } else { 1592 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1593 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1594 1595 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1596 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1597 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1598 /* determine lens of each row */ 1599 for (i=0; i<nrows; i++) { 1600 kstart = ai[irow[i]]; 1601 kend = kstart + a->ilen[irow[i]]; 1602 lens[i] = 0; 1603 for (k=kstart; k<kend; k++) { 1604 if (smap[aj[k]]) { 1605 lens[i]++; 1606 } 1607 } 1608 } 1609 /* Create and fill new matrix */ 1610 if (scall == MAT_REUSE_MATRIX) { 1611 PetscTruth equal; 1612 1613 c = (Mat_SeqAIJ *)((*B)->data); 1614 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1615 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1616 if (!equal) { 1617 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1618 } 1619 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1620 C = *B; 1621 } else { 1622 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1623 } 1624 c = (Mat_SeqAIJ *)(C->data); 1625 for (i=0; i<nrows; i++) { 1626 row = irow[i]; 1627 kstart = ai[row]; 1628 kend = kstart + a->ilen[row]; 1629 mat_i = c->i[i]; 1630 mat_j = c->j + mat_i; 1631 mat_a = c->a + mat_i; 1632 mat_ilen = c->ilen + i; 1633 for (k=kstart; k<kend; k++) { 1634 if ((tcol=smap[a->j[k]])) { 1635 *mat_j++ = tcol - 1; 1636 *mat_a++ = a->a[k]; 1637 (*mat_ilen)++; 1638 1639 } 1640 } 1641 } 1642 /* Free work space */ 1643 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1644 ierr = PetscFree(smap);CHKERRQ(ierr); 1645 ierr = PetscFree(lens);CHKERRQ(ierr); 1646 } 1647 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1648 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1649 1650 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1651 *B = C; 1652 PetscFunctionReturn(0); 1653 } 1654 1655 /* 1656 */ 1657 #undef __FUNCT__ 1658 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1659 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1660 { 1661 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1662 int ierr; 1663 Mat outA; 1664 PetscTruth row_identity,col_identity; 1665 1666 PetscFunctionBegin; 1667 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1668 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1669 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1670 if (!row_identity || !col_identity) { 1671 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1672 } 1673 1674 outA = inA; 1675 inA->factor = FACTOR_LU; 1676 a->row = row; 1677 a->col = col; 1678 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1679 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1680 1681 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1682 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1683 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1684 PetscLogObjectParent(inA,a->icol); 1685 1686 if (!a->solve_work) { /* this matrix may have been factored before */ 1687 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1688 } 1689 1690 if (!a->diag) { 1691 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1692 } 1693 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1694 PetscFunctionReturn(0); 1695 } 1696 1697 #include "petscblaslapack.h" 1698 #undef __FUNCT__ 1699 #define __FUNCT__ "MatScale_SeqAIJ" 1700 int MatScale_SeqAIJ(const PetscScalar *alpha,Mat inA) 1701 { 1702 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1703 int one = 1; 1704 1705 PetscFunctionBegin; 1706 BLscal_(&a->nz,(PetscScalar*)alpha,a->a,&one); 1707 PetscLogFlops(a->nz); 1708 PetscFunctionReturn(0); 1709 } 1710 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1713 int MatGetSubMatrices_SeqAIJ(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1714 { 1715 int ierr,i; 1716 1717 PetscFunctionBegin; 1718 if (scall == MAT_INITIAL_MATRIX) { 1719 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1720 } 1721 1722 for (i=0; i<n; i++) { 1723 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1724 } 1725 PetscFunctionReturn(0); 1726 } 1727 1728 #undef __FUNCT__ 1729 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1730 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1731 { 1732 PetscFunctionBegin; 1733 *bs = 1; 1734 PetscFunctionReturn(0); 1735 } 1736 1737 #undef __FUNCT__ 1738 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1739 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS is[],int ov) 1740 { 1741 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1742 int row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1743 int start,end,*ai,*aj; 1744 PetscBT table; 1745 1746 PetscFunctionBegin; 1747 m = A->m; 1748 ai = a->i; 1749 aj = a->j; 1750 1751 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1752 1753 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1754 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1755 1756 for (i=0; i<is_max; i++) { 1757 /* Initialize the two local arrays */ 1758 isz = 0; 1759 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1760 1761 /* Extract the indices, assume there can be duplicate entries */ 1762 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1763 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1764 1765 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1766 for (j=0; j<n ; ++j){ 1767 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1768 } 1769 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1770 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1771 1772 k = 0; 1773 for (j=0; j<ov; j++){ /* for each overlap */ 1774 n = isz; 1775 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1776 row = nidx[k]; 1777 start = ai[row]; 1778 end = ai[row+1]; 1779 for (l = start; l<end ; l++){ 1780 val = aj[l] ; 1781 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1782 } 1783 } 1784 } 1785 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1786 } 1787 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1788 ierr = PetscFree(nidx);CHKERRQ(ierr); 1789 PetscFunctionReturn(0); 1790 } 1791 1792 /* -------------------------------------------------------------- */ 1793 #undef __FUNCT__ 1794 #define __FUNCT__ "MatPermute_SeqAIJ" 1795 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1796 { 1797 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1798 PetscScalar *vwork; 1799 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1800 int *row,*col,*cnew,j,*lens; 1801 IS icolp,irowp; 1802 1803 PetscFunctionBegin; 1804 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1805 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1806 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1807 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1808 1809 /* determine lengths of permuted rows */ 1810 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1811 for (i=0; i<m; i++) { 1812 lens[row[i]] = a->i[i+1] - a->i[i]; 1813 } 1814 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1815 ierr = PetscFree(lens);CHKERRQ(ierr); 1816 1817 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1818 for (i=0; i<m; i++) { 1819 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1820 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1821 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1822 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1823 } 1824 ierr = PetscFree(cnew);CHKERRQ(ierr); 1825 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1826 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1827 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1828 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1829 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1830 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1831 PetscFunctionReturn(0); 1832 } 1833 1834 #undef __FUNCT__ 1835 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1836 int MatPrintHelp_SeqAIJ(Mat A) 1837 { 1838 static PetscTruth called = PETSC_FALSE; 1839 MPI_Comm comm = A->comm; 1840 int ierr; 1841 1842 PetscFunctionBegin; 1843 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1844 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1845 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1846 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1847 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1848 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1849 #if defined(PETSC_HAVE_ESSL) 1850 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1851 #endif 1852 #if defined(PETSC_HAVE_LUSOL) 1853 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1854 #endif 1855 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1856 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1857 #endif 1858 PetscFunctionReturn(0); 1859 } 1860 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1861 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1862 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatFactorInfo*,IS,IS,Mat*); 1863 #undef __FUNCT__ 1864 #define __FUNCT__ "MatCopy_SeqAIJ" 1865 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1866 { 1867 int ierr; 1868 1869 PetscFunctionBegin; 1870 /* If the two matrices have the same copy implementation, use fast copy. */ 1871 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1872 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1873 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1874 1875 if (a->i[A->m] != b->i[B->m]) { 1876 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1877 } 1878 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1879 } else { 1880 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1881 } 1882 PetscFunctionReturn(0); 1883 } 1884 1885 #undef __FUNCT__ 1886 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1887 int MatSetUpPreallocation_SeqAIJ(Mat A) 1888 { 1889 int ierr; 1890 1891 PetscFunctionBegin; 1892 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1893 PetscFunctionReturn(0); 1894 } 1895 1896 #undef __FUNCT__ 1897 #define __FUNCT__ "MatGetArray_SeqAIJ" 1898 int MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1899 { 1900 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1901 PetscFunctionBegin; 1902 *array = a->a; 1903 PetscFunctionReturn(0); 1904 } 1905 1906 #undef __FUNCT__ 1907 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1908 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1909 { 1910 PetscFunctionBegin; 1911 PetscFunctionReturn(0); 1912 } 1913 1914 #undef __FUNCT__ 1915 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1916 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1917 { 1918 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1919 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1920 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1921 PetscScalar *vscale_array; 1922 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1923 Vec w1,w2,w3; 1924 void *fctx = coloring->fctx; 1925 PetscTruth flg; 1926 1927 PetscFunctionBegin; 1928 if (!coloring->w1) { 1929 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1930 PetscLogObjectParent(coloring,coloring->w1); 1931 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1932 PetscLogObjectParent(coloring,coloring->w2); 1933 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1934 PetscLogObjectParent(coloring,coloring->w3); 1935 } 1936 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1937 1938 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1939 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1940 if (flg) { 1941 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1942 } else { 1943 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1944 } 1945 1946 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1947 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1948 1949 /* 1950 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1951 coloring->F for the coarser grids from the finest 1952 */ 1953 if (coloring->F) { 1954 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1955 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1956 if (m1 != m2) { 1957 coloring->F = 0; 1958 } 1959 } 1960 1961 if (coloring->F) { 1962 w1 = coloring->F; 1963 coloring->F = 0; 1964 } else { 1965 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1966 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1967 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1968 } 1969 1970 /* 1971 Compute all the scale factors and share with other processors 1972 */ 1973 ierr = VecGetArrayFast(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1974 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1975 for (k=0; k<coloring->ncolors; k++) { 1976 /* 1977 Loop over each column associated with color adding the 1978 perturbation to the vector w3. 1979 */ 1980 for (l=0; l<coloring->ncolumns[k]; l++) { 1981 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1982 dx = xx[col]; 1983 if (dx == 0.0) dx = 1.0; 1984 #if !defined(PETSC_USE_COMPLEX) 1985 if (dx < umin && dx >= 0.0) dx = umin; 1986 else if (dx < 0.0 && dx > -umin) dx = -umin; 1987 #else 1988 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1989 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1990 #endif 1991 dx *= epsilon; 1992 vscale_array[col] = 1.0/dx; 1993 } 1994 } 1995 vscale_array = vscale_array + start;ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1996 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1997 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1998 1999 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2000 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2001 2002 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2003 else vscaleforrow = coloring->columnsforrow; 2004 2005 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2006 /* 2007 Loop over each color 2008 */ 2009 for (k=0; k<coloring->ncolors; k++) { 2010 coloring->currentcolor = k; 2011 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2012 ierr = VecGetArrayFast(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2013 /* 2014 Loop over each column associated with color adding the 2015 perturbation to the vector w3. 2016 */ 2017 for (l=0; l<coloring->ncolumns[k]; l++) { 2018 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2019 dx = xx[col]; 2020 if (dx == 0.0) dx = 1.0; 2021 #if !defined(PETSC_USE_COMPLEX) 2022 if (dx < umin && dx >= 0.0) dx = umin; 2023 else if (dx < 0.0 && dx > -umin) dx = -umin; 2024 #else 2025 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2026 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2027 #endif 2028 dx *= epsilon; 2029 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2030 w3_array[col] += dx; 2031 } 2032 w3_array = w3_array + start; ierr = VecRestoreArrayFast(w3,&w3_array);CHKERRQ(ierr); 2033 2034 /* 2035 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2036 */ 2037 2038 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2039 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2040 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2041 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2042 2043 /* 2044 Loop over rows of vector, putting results into Jacobian matrix 2045 */ 2046 ierr = VecGetArrayFast(w2,&y);CHKERRQ(ierr); 2047 for (l=0; l<coloring->nrows[k]; l++) { 2048 row = coloring->rows[k][l]; 2049 col = coloring->columnsforrow[k][l]; 2050 y[row] *= vscale_array[vscaleforrow[k][l]]; 2051 srow = row + start; 2052 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2053 } 2054 ierr = VecRestoreArrayFast(w2,&y);CHKERRQ(ierr); 2055 } 2056 coloring->currentcolor = k; 2057 ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2058 xx = xx + start; ierr = VecRestoreArrayFast(x1,&xx);CHKERRQ(ierr); 2059 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2060 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2061 PetscFunctionReturn(0); 2062 } 2063 2064 #include "petscblaslapack.h" 2065 #undef __FUNCT__ 2066 #define __FUNCT__ "MatAXPY_SeqAIJ" 2067 int MatAXPY_SeqAIJ(const PetscScalar *a,Mat X,Mat Y,MatStructure str) 2068 { 2069 int ierr,one=1,i; 2070 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2071 2072 PetscFunctionBegin; 2073 if (str == SAME_NONZERO_PATTERN) { 2074 BLaxpy_(&x->nz,(PetscScalar*)a,x->a,&one,y->a,&one); 2075 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2076 if (y->xtoy && y->XtoY != X) { 2077 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2078 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2079 } 2080 if (!y->xtoy) { /* get xtoy */ 2081 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2082 y->XtoY = X; 2083 } 2084 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2085 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2086 } else { 2087 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2088 } 2089 PetscFunctionReturn(0); 2090 } 2091 2092 /* -------------------------------------------------------------------*/ 2093 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2094 MatGetRow_SeqAIJ, 2095 MatRestoreRow_SeqAIJ, 2096 MatMult_SeqAIJ, 2097 MatMultAdd_SeqAIJ, 2098 MatMultTranspose_SeqAIJ, 2099 MatMultTransposeAdd_SeqAIJ, 2100 MatSolve_SeqAIJ, 2101 MatSolveAdd_SeqAIJ, 2102 MatSolveTranspose_SeqAIJ, 2103 MatSolveTransposeAdd_SeqAIJ, 2104 MatLUFactor_SeqAIJ, 2105 0, 2106 MatRelax_SeqAIJ, 2107 MatTranspose_SeqAIJ, 2108 MatGetInfo_SeqAIJ, 2109 MatEqual_SeqAIJ, 2110 MatGetDiagonal_SeqAIJ, 2111 MatDiagonalScale_SeqAIJ, 2112 MatNorm_SeqAIJ, 2113 0, 2114 MatAssemblyEnd_SeqAIJ, 2115 MatCompress_SeqAIJ, 2116 MatSetOption_SeqAIJ, 2117 MatZeroEntries_SeqAIJ, 2118 MatZeroRows_SeqAIJ, 2119 MatLUFactorSymbolic_SeqAIJ, 2120 MatLUFactorNumeric_SeqAIJ, 2121 MatCholeskyFactorSymbolic_SeqAIJ, 2122 MatCholeskyFactorNumeric_SeqAIJ, 2123 MatSetUpPreallocation_SeqAIJ, 2124 MatILUFactorSymbolic_SeqAIJ, 2125 MatICCFactorSymbolic_SeqAIJ, 2126 MatGetArray_SeqAIJ, 2127 MatRestoreArray_SeqAIJ, 2128 MatDuplicate_SeqAIJ, 2129 0, 2130 0, 2131 MatILUFactor_SeqAIJ, 2132 0, 2133 MatAXPY_SeqAIJ, 2134 MatGetSubMatrices_SeqAIJ, 2135 MatIncreaseOverlap_SeqAIJ, 2136 MatGetValues_SeqAIJ, 2137 MatCopy_SeqAIJ, 2138 MatPrintHelp_SeqAIJ, 2139 MatScale_SeqAIJ, 2140 0, 2141 0, 2142 MatILUDTFactor_SeqAIJ, 2143 MatGetBlockSize_SeqAIJ, 2144 MatGetRowIJ_SeqAIJ, 2145 MatRestoreRowIJ_SeqAIJ, 2146 MatGetColumnIJ_SeqAIJ, 2147 MatRestoreColumnIJ_SeqAIJ, 2148 MatFDColoringCreate_SeqAIJ, 2149 0, 2150 0, 2151 MatPermute_SeqAIJ, 2152 0, 2153 0, 2154 MatDestroy_SeqAIJ, 2155 MatView_SeqAIJ, 2156 MatGetPetscMaps_Petsc, 2157 0, 2158 0, 2159 0, 2160 0, 2161 0, 2162 0, 2163 0, 2164 0, 2165 MatSetColoring_SeqAIJ, 2166 MatSetValuesAdic_SeqAIJ, 2167 MatSetValuesAdifor_SeqAIJ, 2168 MatFDColoringApply_SeqAIJ}; 2169 2170 EXTERN_C_BEGIN 2171 #undef __FUNCT__ 2172 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2173 2174 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2175 { 2176 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2177 int i,nz,n; 2178 2179 PetscFunctionBegin; 2180 2181 nz = aij->maxnz; 2182 n = mat->n; 2183 for (i=0; i<nz; i++) { 2184 aij->j[i] = indices[i]; 2185 } 2186 aij->nz = nz; 2187 for (i=0; i<n; i++) { 2188 aij->ilen[i] = aij->imax[i]; 2189 } 2190 2191 PetscFunctionReturn(0); 2192 } 2193 EXTERN_C_END 2194 2195 #undef __FUNCT__ 2196 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2197 /*@ 2198 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2199 in the matrix. 2200 2201 Input Parameters: 2202 + mat - the SeqAIJ matrix 2203 - indices - the column indices 2204 2205 Level: advanced 2206 2207 Notes: 2208 This can be called if you have precomputed the nonzero structure of the 2209 matrix and want to provide it to the matrix object to improve the performance 2210 of the MatSetValues() operation. 2211 2212 You MUST have set the correct numbers of nonzeros per row in the call to 2213 MatCreateSeqAIJ(). 2214 2215 MUST be called before any calls to MatSetValues(); 2216 2217 The indices should start with zero, not one. 2218 2219 @*/ 2220 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2221 { 2222 int ierr,(*f)(Mat,int *); 2223 2224 PetscFunctionBegin; 2225 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2226 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2227 if (f) { 2228 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2229 } else { 2230 SETERRQ(1,"Wrong type of matrix to set column indices"); 2231 } 2232 PetscFunctionReturn(0); 2233 } 2234 2235 /* ----------------------------------------------------------------------------------------*/ 2236 2237 EXTERN_C_BEGIN 2238 #undef __FUNCT__ 2239 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2240 int MatStoreValues_SeqAIJ(Mat mat) 2241 { 2242 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2243 size_t nz = aij->i[mat->m],ierr; 2244 2245 PetscFunctionBegin; 2246 if (aij->nonew != 1) { 2247 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2248 } 2249 2250 /* allocate space for values if not already there */ 2251 if (!aij->saved_values) { 2252 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2253 } 2254 2255 /* copy values over */ 2256 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2257 PetscFunctionReturn(0); 2258 } 2259 EXTERN_C_END 2260 2261 #undef __FUNCT__ 2262 #define __FUNCT__ "MatStoreValues" 2263 /*@ 2264 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2265 example, reuse of the linear part of a Jacobian, while recomputing the 2266 nonlinear portion. 2267 2268 Collect on Mat 2269 2270 Input Parameters: 2271 . mat - the matrix (currently on AIJ matrices support this option) 2272 2273 Level: advanced 2274 2275 Common Usage, with SNESSolve(): 2276 $ Create Jacobian matrix 2277 $ Set linear terms into matrix 2278 $ Apply boundary conditions to matrix, at this time matrix must have 2279 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2280 $ boundary conditions again will not change the nonzero structure 2281 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2282 $ ierr = MatStoreValues(mat); 2283 $ Call SNESSetJacobian() with matrix 2284 $ In your Jacobian routine 2285 $ ierr = MatRetrieveValues(mat); 2286 $ Set nonlinear terms in matrix 2287 2288 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2289 $ // build linear portion of Jacobian 2290 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2291 $ ierr = MatStoreValues(mat); 2292 $ loop over nonlinear iterations 2293 $ ierr = MatRetrieveValues(mat); 2294 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2295 $ // call MatAssemblyBegin/End() on matrix 2296 $ Solve linear system with Jacobian 2297 $ endloop 2298 2299 Notes: 2300 Matrix must already be assemblied before calling this routine 2301 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2302 calling this routine. 2303 2304 .seealso: MatRetrieveValues() 2305 2306 @*/ 2307 int MatStoreValues(Mat mat) 2308 { 2309 int ierr,(*f)(Mat); 2310 2311 PetscFunctionBegin; 2312 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2313 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2314 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2315 2316 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2317 if (f) { 2318 ierr = (*f)(mat);CHKERRQ(ierr); 2319 } else { 2320 SETERRQ(1,"Wrong type of matrix to store values"); 2321 } 2322 PetscFunctionReturn(0); 2323 } 2324 2325 EXTERN_C_BEGIN 2326 #undef __FUNCT__ 2327 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2328 int MatRetrieveValues_SeqAIJ(Mat mat) 2329 { 2330 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2331 int nz = aij->i[mat->m],ierr; 2332 2333 PetscFunctionBegin; 2334 if (aij->nonew != 1) { 2335 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2336 } 2337 if (!aij->saved_values) { 2338 SETERRQ(1,"Must call MatStoreValues(A);first"); 2339 } 2340 2341 /* copy values over */ 2342 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2343 PetscFunctionReturn(0); 2344 } 2345 EXTERN_C_END 2346 2347 #undef __FUNCT__ 2348 #define __FUNCT__ "MatRetrieveValues" 2349 /*@ 2350 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2351 example, reuse of the linear part of a Jacobian, while recomputing the 2352 nonlinear portion. 2353 2354 Collect on Mat 2355 2356 Input Parameters: 2357 . mat - the matrix (currently on AIJ matrices support this option) 2358 2359 Level: advanced 2360 2361 .seealso: MatStoreValues() 2362 2363 @*/ 2364 int MatRetrieveValues(Mat mat) 2365 { 2366 int ierr,(*f)(Mat); 2367 2368 PetscFunctionBegin; 2369 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2370 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2371 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2372 2373 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2374 if (f) { 2375 ierr = (*f)(mat);CHKERRQ(ierr); 2376 } else { 2377 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2378 } 2379 PetscFunctionReturn(0); 2380 } 2381 2382 /* 2383 This allows SeqAIJ matrices to be passed to the matlab engine 2384 */ 2385 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2386 #include "engine.h" /* Matlab include file */ 2387 #include "mex.h" /* Matlab include file */ 2388 EXTERN_C_BEGIN 2389 #undef __FUNCT__ 2390 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2391 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 2392 { 2393 int ierr,i,*ai,*aj; 2394 Mat B = (Mat)obj; 2395 mxArray *mat; 2396 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2397 2398 PetscFunctionBegin; 2399 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2400 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2401 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2402 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2403 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2404 2405 /* Matlab indices start at 0 for sparse (what a surprise) */ 2406 2407 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2408 mxSetName(mat,obj->name); 2409 engPutArray((Engine *)mengine,mat); 2410 PetscFunctionReturn(0); 2411 } 2412 EXTERN_C_END 2413 2414 EXTERN_C_BEGIN 2415 #undef __FUNCT__ 2416 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2417 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 2418 { 2419 int ierr,ii; 2420 Mat mat = (Mat)obj; 2421 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2422 mxArray *mmat; 2423 2424 PetscFunctionBegin; 2425 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2426 2427 mmat = engGetArray((Engine *)mengine,obj->name); 2428 2429 aij->nz = (mxGetJc(mmat))[mat->m]; 2430 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2431 aij->j = (int*)(aij->a + aij->nz); 2432 aij->i = aij->j + aij->nz; 2433 aij->singlemalloc = PETSC_TRUE; 2434 aij->freedata = PETSC_TRUE; 2435 2436 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2437 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2438 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2439 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2440 2441 for (ii=0; ii<mat->m; ii++) { 2442 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2443 } 2444 2445 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2446 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2447 2448 PetscFunctionReturn(0); 2449 } 2450 EXTERN_C_END 2451 #endif 2452 2453 /* --------------------------------------------------------------------------------*/ 2454 #undef __FUNCT__ 2455 #define __FUNCT__ "MatCreateSeqAIJ" 2456 /*@C 2457 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2458 (the default parallel PETSc format). For good matrix assembly performance 2459 the user should preallocate the matrix storage by setting the parameter nz 2460 (or the array nnz). By setting these parameters accurately, performance 2461 during matrix assembly can be increased by more than a factor of 50. 2462 2463 Collective on MPI_Comm 2464 2465 Input Parameters: 2466 + comm - MPI communicator, set to PETSC_COMM_SELF 2467 . m - number of rows 2468 . n - number of columns 2469 . nz - number of nonzeros per row (same for all rows) 2470 - nnz - array containing the number of nonzeros in the various rows 2471 (possibly different for each row) or PETSC_NULL 2472 2473 Output Parameter: 2474 . A - the matrix 2475 2476 Notes: 2477 The AIJ format (also called the Yale sparse matrix format or 2478 compressed row storage), is fully compatible with standard Fortran 77 2479 storage. That is, the stored row and column indices can begin at 2480 either one (as in Fortran) or zero. See the users' manual for details. 2481 2482 Specify the preallocated storage with either nz or nnz (not both). 2483 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2484 allocation. For large problems you MUST preallocate memory or you 2485 will get TERRIBLE performance, see the users' manual chapter on matrices. 2486 2487 By default, this format uses inodes (identical nodes) when possible, to 2488 improve numerical efficiency of matrix-vector products and solves. We 2489 search for consecutive rows with the same nonzero structure, thereby 2490 reusing matrix information to achieve increased efficiency. 2491 2492 Options Database Keys: 2493 + -mat_aij_no_inode - Do not use inodes 2494 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2495 - -mat_aij_oneindex - Internally use indexing starting at 1 2496 rather than 0. Note that when calling MatSetValues(), 2497 the user still MUST index entries starting at 0! 2498 2499 Level: intermediate 2500 2501 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2502 2503 @*/ 2504 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,const int nnz[],Mat *A) 2505 { 2506 int ierr; 2507 2508 PetscFunctionBegin; 2509 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2510 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2511 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2512 PetscFunctionReturn(0); 2513 } 2514 2515 #define SKIP_ALLOCATION -4 2516 2517 #undef __FUNCT__ 2518 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2519 /*@C 2520 MatSeqAIJSetPreallocation - For good matrix assembly performance 2521 the user should preallocate the matrix storage by setting the parameter nz 2522 (or the array nnz). By setting these parameters accurately, performance 2523 during matrix assembly can be increased by more than a factor of 50. 2524 2525 Collective on MPI_Comm 2526 2527 Input Parameters: 2528 + comm - MPI communicator, set to PETSC_COMM_SELF 2529 . m - number of rows 2530 . n - number of columns 2531 . nz - number of nonzeros per row (same for all rows) 2532 - nnz - array containing the number of nonzeros in the various rows 2533 (possibly different for each row) or PETSC_NULL 2534 2535 Output Parameter: 2536 . A - the matrix 2537 2538 Notes: 2539 The AIJ format (also called the Yale sparse matrix format or 2540 compressed row storage), is fully compatible with standard Fortran 77 2541 storage. That is, the stored row and column indices can begin at 2542 either one (as in Fortran) or zero. See the users' manual for details. 2543 2544 Specify the preallocated storage with either nz or nnz (not both). 2545 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2546 allocation. For large problems you MUST preallocate memory or you 2547 will get TERRIBLE performance, see the users' manual chapter on matrices. 2548 2549 By default, this format uses inodes (identical nodes) when possible, to 2550 improve numerical efficiency of matrix-vector products and solves. We 2551 search for consecutive rows with the same nonzero structure, thereby 2552 reusing matrix information to achieve increased efficiency. 2553 2554 Options Database Keys: 2555 + -mat_aij_no_inode - Do not use inodes 2556 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2557 - -mat_aij_oneindex - Internally use indexing starting at 1 2558 rather than 0. Note that when calling MatSetValues(), 2559 the user still MUST index entries starting at 0! 2560 2561 Level: intermediate 2562 2563 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2564 2565 @*/ 2566 int MatSeqAIJSetPreallocation(Mat B,int nz,const int nnz[]) 2567 { 2568 int ierr,(*f)(Mat,int,const int[]); 2569 2570 PetscFunctionBegin; 2571 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2572 if (f) { 2573 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2574 } 2575 PetscFunctionReturn(0); 2576 } 2577 2578 EXTERN_C_BEGIN 2579 #undef __FUNCT__ 2580 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2581 int MatSeqAIJSetPreallocation_SeqAIJ(Mat B,int nz,int *nnz) 2582 { 2583 Mat_SeqAIJ *b; 2584 size_t len = 0; 2585 PetscTruth skipallocation = PETSC_FALSE; 2586 int i,ierr; 2587 2588 PetscFunctionBegin; 2589 2590 if (nz == SKIP_ALLOCATION) { 2591 skipallocation = PETSC_TRUE; 2592 nz = 0; 2593 } 2594 2595 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2596 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2597 if (nnz) { 2598 for (i=0; i<B->m; i++) { 2599 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2600 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); 2601 } 2602 } 2603 2604 B->preallocated = PETSC_TRUE; 2605 b = (Mat_SeqAIJ*)B->data; 2606 2607 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2608 if (!nnz) { 2609 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2610 else if (nz <= 0) nz = 1; 2611 for (i=0; i<B->m; i++) b->imax[i] = nz; 2612 nz = nz*B->m; 2613 } else { 2614 nz = 0; 2615 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2616 } 2617 2618 if (!skipallocation) { 2619 /* allocate the matrix space */ 2620 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2621 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2622 b->j = (int*)(b->a + nz); 2623 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2624 b->i = b->j + nz; 2625 b->i[0] = 0; 2626 for (i=1; i<B->m+1; i++) { 2627 b->i[i] = b->i[i-1] + b->imax[i-1]; 2628 } 2629 b->singlemalloc = PETSC_TRUE; 2630 b->freedata = PETSC_TRUE; 2631 } else { 2632 b->freedata = PETSC_FALSE; 2633 } 2634 2635 /* b->ilen will count nonzeros in each row so far. */ 2636 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2637 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2638 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2639 2640 b->nz = 0; 2641 b->maxnz = nz; 2642 B->info.nz_unneeded = (double)b->maxnz; 2643 PetscFunctionReturn(0); 2644 } 2645 EXTERN_C_END 2646 2647 EXTERN int RegisterApplyPtAPRoutines_Private(Mat); 2648 2649 EXTERN_C_BEGIN 2650 EXTERN int MatConvert_SeqAIJ_SeqSBAIJ(Mat,MatType,Mat*); 2651 EXTERN int MatConvert_SeqAIJ_SeqBAIJ(Mat,MatType,Mat*); 2652 EXTERN int MatReorderForNonzeroDiagonal_SeqAIJ(Mat,PetscReal,IS,IS); 2653 EXTERN int MatAdjustForInodes_SeqAIJ(Mat,IS*,IS*); 2654 EXTERN int MatSeqAIJGetInodeSizes_SeqAIJ(Mat,int*,int*[],int*); 2655 EXTERN_C_END 2656 2657 EXTERN_C_BEGIN 2658 #undef __FUNCT__ 2659 #define __FUNCT__ "MatCreate_SeqAIJ" 2660 int MatCreate_SeqAIJ(Mat B) 2661 { 2662 Mat_SeqAIJ *b; 2663 int ierr,size; 2664 PetscTruth flg; 2665 2666 PetscFunctionBegin; 2667 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2668 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2669 2670 B->m = B->M = PetscMax(B->m,B->M); 2671 B->n = B->N = PetscMax(B->n,B->N); 2672 2673 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2674 B->data = (void*)b; 2675 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2676 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2677 B->factor = 0; 2678 B->lupivotthreshold = 1.0; 2679 B->mapping = 0; 2680 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2681 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2682 b->row = 0; 2683 b->col = 0; 2684 b->icol = 0; 2685 b->reallocs = 0; 2686 2687 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2688 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2689 2690 b->sorted = PETSC_FALSE; 2691 b->ignorezeroentries = PETSC_FALSE; 2692 b->roworiented = PETSC_TRUE; 2693 b->nonew = 0; 2694 b->diag = 0; 2695 b->solve_work = 0; 2696 B->spptr = 0; 2697 b->inode.use = PETSC_TRUE; 2698 b->inode.node_count = 0; 2699 b->inode.size = 0; 2700 b->inode.limit = 5; 2701 b->inode.max_limit = 5; 2702 b->saved_values = 0; 2703 b->idiag = 0; 2704 b->ssor = 0; 2705 b->keepzeroedrows = PETSC_FALSE; 2706 b->xtoy = 0; 2707 b->XtoY = 0; 2708 2709 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2710 2711 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2712 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2713 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2714 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2715 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2716 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2717 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2718 if (flg){ 2719 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2720 } 2721 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2722 "MatSeqAIJSetColumnIndices_SeqAIJ", 2723 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2724 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2725 "MatStoreValues_SeqAIJ", 2726 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2727 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2728 "MatRetrieveValues_SeqAIJ", 2729 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2730 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2731 "MatConvert_SeqAIJ_SeqSBAIJ", 2732 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2733 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2734 "MatConvert_SeqAIJ_SeqBAIJ", 2735 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2736 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsSymmetric_C", 2737 "MatIsSymmetric_SeqAIJ", 2738 MatIsSymmetric_SeqAIJ);CHKERRQ(ierr); 2739 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2740 "MatSeqAIJSetPreallocation_SeqAIJ", 2741 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2742 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2743 "MatReorderForNonzeroDiagonal_SeqAIJ", 2744 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2745 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2746 "MatAdjustForInodes_SeqAIJ", 2747 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2748 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2749 "MatSeqAIJGetInodeSizes_SeqAIJ", 2750 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2751 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2752 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2753 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2754 #endif 2755 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2756 PetscFunctionReturn(0); 2757 } 2758 EXTERN_C_END 2759 2760 #undef __FUNCT__ 2761 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2762 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2763 { 2764 Mat C; 2765 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2766 int i,m = A->m,ierr; 2767 size_t len; 2768 2769 PetscFunctionBegin; 2770 *B = 0; 2771 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2772 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2773 c = (Mat_SeqAIJ*)C->data; 2774 2775 C->factor = A->factor; 2776 c->row = 0; 2777 c->col = 0; 2778 c->icol = 0; 2779 c->keepzeroedrows = a->keepzeroedrows; 2780 C->assembled = PETSC_TRUE; 2781 2782 C->M = A->m; 2783 C->N = A->n; 2784 2785 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2786 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2787 for (i=0; i<m; i++) { 2788 c->imax[i] = a->imax[i]; 2789 c->ilen[i] = a->ilen[i]; 2790 } 2791 2792 /* allocate the matrix space */ 2793 c->singlemalloc = PETSC_TRUE; 2794 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2795 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2796 c->j = (int*)(c->a + a->i[m] ); 2797 c->i = c->j + a->i[m]; 2798 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2799 if (m > 0) { 2800 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(int));CHKERRQ(ierr); 2801 if (cpvalues == MAT_COPY_VALUES) { 2802 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2803 } else { 2804 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2805 } 2806 } 2807 2808 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2809 c->sorted = a->sorted; 2810 c->roworiented = a->roworiented; 2811 c->nonew = a->nonew; 2812 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2813 c->saved_values = 0; 2814 c->idiag = 0; 2815 c->ssor = 0; 2816 c->ignorezeroentries = a->ignorezeroentries; 2817 c->freedata = PETSC_TRUE; 2818 2819 if (a->diag) { 2820 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2821 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2822 for (i=0; i<m; i++) { 2823 c->diag[i] = a->diag[i]; 2824 } 2825 } else c->diag = 0; 2826 c->inode.use = a->inode.use; 2827 c->inode.limit = a->inode.limit; 2828 c->inode.max_limit = a->inode.max_limit; 2829 if (a->inode.size){ 2830 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2831 c->inode.node_count = a->inode.node_count; 2832 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2833 } else { 2834 c->inode.size = 0; 2835 c->inode.node_count = 0; 2836 } 2837 c->nz = a->nz; 2838 c->maxnz = a->maxnz; 2839 c->solve_work = 0; 2840 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2841 C->preallocated = PETSC_TRUE; 2842 2843 *B = C; 2844 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2845 PetscFunctionReturn(0); 2846 } 2847 2848 EXTERN_C_BEGIN 2849 #undef __FUNCT__ 2850 #define __FUNCT__ "MatLoad_SeqAIJ" 2851 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2852 { 2853 Mat_SeqAIJ *a; 2854 Mat B; 2855 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N; 2856 MPI_Comm comm; 2857 #if defined(PETSC_HAVE_SUPERLUDIST) || defined(PETSC_HAVE_UMFPACK) || defined(PETSC_HAVE_MUMPS) 2858 PetscTruth flag; 2859 #endif 2860 2861 PetscFunctionBegin; 2862 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2863 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2864 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2865 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2866 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2867 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2868 M = header[1]; N = header[2]; nz = header[3]; 2869 2870 if (nz < 0) { 2871 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2872 } 2873 2874 /* read in row lengths */ 2875 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2876 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2877 2878 /* create our matrix */ 2879 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2880 ierr = MatSetType(B,type);CHKERRQ(ierr); 2881 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2882 a = (Mat_SeqAIJ*)B->data; 2883 2884 /* read in column indices and adjust for Fortran indexing*/ 2885 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2886 2887 /* read in nonzero values */ 2888 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2889 2890 /* set matrix "i" values */ 2891 a->i[0] = 0; 2892 for (i=1; i<= M; i++) { 2893 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2894 a->ilen[i-1] = rowlengths[i-1]; 2895 } 2896 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2897 2898 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2899 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2900 #if defined(PETSC_HAVE_SUPERLUDIST) 2901 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_superlu_dist",&flag);CHKERRQ(ierr); 2902 if (flag) { ierr = MatUseSuperLU_DIST_MPIAIJ(B);CHKERRQ(ierr); } 2903 #endif 2904 #if defined(PETSC_HAVE_UMFPACK) 2905 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_umfpack",&flag);CHKERRQ(ierr); 2906 if (flag) { ierr = MatUseUMFPACK_SeqAIJ(B);CHKERRQ(ierr); } 2907 #endif 2908 #if defined(PETSC_HAVE_MUMPS) 2909 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_mumps",&flag);CHKERRQ(ierr); 2910 if (flag) { ierr = MatUseMUMPS_MPIAIJ(B);CHKERRQ(ierr); } 2911 #endif 2912 *A = B; 2913 PetscFunctionReturn(0); 2914 } 2915 EXTERN_C_END 2916 2917 #undef __FUNCT__ 2918 #define __FUNCT__ "MatEqual_SeqAIJ" 2919 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2920 { 2921 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2922 int ierr; 2923 2924 PetscFunctionBegin; 2925 /* If the matrix dimensions are not equal,or no of nonzeros */ 2926 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2927 *flg = PETSC_FALSE; 2928 PetscFunctionReturn(0); 2929 } 2930 2931 /* if the a->i are the same */ 2932 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2933 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2934 2935 /* if a->j are the same */ 2936 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2937 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2938 2939 /* if a->a are the same */ 2940 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2941 2942 PetscFunctionReturn(0); 2943 2944 } 2945 2946 #undef __FUNCT__ 2947 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2948 /*@C 2949 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2950 provided by the user. 2951 2952 Coolective on MPI_Comm 2953 2954 Input Parameters: 2955 + comm - must be an MPI communicator of size 1 2956 . m - number of rows 2957 . n - number of columns 2958 . i - row indices 2959 . j - column indices 2960 - a - matrix values 2961 2962 Output Parameter: 2963 . mat - the matrix 2964 2965 Level: intermediate 2966 2967 Notes: 2968 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2969 once the matrix is destroyed 2970 2971 You cannot set new nonzero locations into this matrix, that will generate an error. 2972 2973 The i and j indices are 0 based 2974 2975 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2976 2977 @*/ 2978 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2979 { 2980 int ierr,ii; 2981 Mat_SeqAIJ *aij; 2982 2983 PetscFunctionBegin; 2984 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2985 aij = (Mat_SeqAIJ*)(*mat)->data; 2986 2987 if (i[0] != 0) { 2988 SETERRQ(1,"i (row indices) must start with 0"); 2989 } 2990 aij->i = i; 2991 aij->j = j; 2992 aij->a = a; 2993 aij->singlemalloc = PETSC_FALSE; 2994 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2995 aij->freedata = PETSC_FALSE; 2996 2997 for (ii=0; ii<m; ii++) { 2998 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2999 #if defined(PETSC_USE_BOPT_g) 3000 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]); 3001 #endif 3002 } 3003 #if defined(PETSC_USE_BOPT_g) 3004 for (ii=0; ii<aij->i[m]; ii++) { 3005 if (j[ii] < 0) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 3006 if (j[ii] > n - 1) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 3007 } 3008 #endif 3009 3010 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3011 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3012 PetscFunctionReturn(0); 3013 } 3014 3015 #undef __FUNCT__ 3016 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3017 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3018 { 3019 int ierr; 3020 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3021 3022 PetscFunctionBegin; 3023 if (coloring->ctype == IS_COLORING_LOCAL) { 3024 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3025 a->coloring = coloring; 3026 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3027 int i,*larray; 3028 ISColoring ocoloring; 3029 ISColoringValue *colors; 3030 3031 /* set coloring for diagonal portion */ 3032 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 3033 for (i=0; i<A->n; i++) { 3034 larray[i] = i; 3035 } 3036 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3037 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3038 for (i=0; i<A->n; i++) { 3039 colors[i] = coloring->colors[larray[i]]; 3040 } 3041 ierr = PetscFree(larray);CHKERRQ(ierr); 3042 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3043 a->coloring = ocoloring; 3044 } 3045 PetscFunctionReturn(0); 3046 } 3047 3048 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3049 EXTERN_C_BEGIN 3050 #include "adic/ad_utils.h" 3051 EXTERN_C_END 3052 3053 #undef __FUNCT__ 3054 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3055 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3056 { 3057 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3058 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3059 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3060 ISColoringValue *color; 3061 3062 PetscFunctionBegin; 3063 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3064 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3065 color = a->coloring->colors; 3066 /* loop over rows */ 3067 for (i=0; i<m; i++) { 3068 nz = ii[i+1] - ii[i]; 3069 /* loop over columns putting computed value into matrix */ 3070 for (j=0; j<nz; j++) { 3071 *v++ = values[color[*jj++]]; 3072 } 3073 values += nlen; /* jump to next row of derivatives */ 3074 } 3075 PetscFunctionReturn(0); 3076 } 3077 3078 #else 3079 3080 #undef __FUNCT__ 3081 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3082 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3083 { 3084 PetscFunctionBegin; 3085 SETERRQ(1,"PETSc installed without ADIC"); 3086 } 3087 3088 #endif 3089 3090 #undef __FUNCT__ 3091 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3092 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3093 { 3094 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3095 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3096 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3097 ISColoringValue *color; 3098 3099 PetscFunctionBegin; 3100 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3101 color = a->coloring->colors; 3102 /* loop over rows */ 3103 for (i=0; i<m; i++) { 3104 nz = ii[i+1] - ii[i]; 3105 /* loop over columns putting computed value into matrix */ 3106 for (j=0; j<nz; j++) { 3107 *v++ = values[color[*jj++]]; 3108 } 3109 values += nl; /* jump to next row of derivatives */ 3110 } 3111 PetscFunctionReturn(0); 3112 } 3113 3114