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