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