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