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