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