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