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