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