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