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