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