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