1 #define PETSCMAT_DLL 2 3 4 /* 5 Defines the basic matrix operations for the AIJ (compressed row) 6 matrix storage format. 7 */ 8 9 10 #include "../src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 11 #include "petscblaslapack.h" 12 #include "petscbt.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 16 PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 17 { 18 PetscErrorCode ierr; 19 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 20 PetscInt i,*diag, m = Y->rmap->n; 21 MatScalar *aa = aij->a; 22 PetscScalar *v; 23 PetscTruth missing; 24 25 PetscFunctionBegin; 26 if (Y->assembled) { 27 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,PETSC_NULL);CHKERRQ(ierr); 28 if (!missing) { 29 diag = aij->diag; 30 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 31 if (is == INSERT_VALUES) { 32 for (i=0; i<m; i++) { 33 aa[diag[i]] = v[i]; 34 } 35 } else { 36 for (i=0; i<m; i++) { 37 aa[diag[i]] += v[i]; 38 } 39 } 40 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 } 44 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 50 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 51 { 52 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 53 PetscErrorCode ierr; 54 PetscInt i,ishift; 55 56 PetscFunctionBegin; 57 *m = A->rmap->n; 58 if (!ia) PetscFunctionReturn(0); 59 ishift = 0; 60 if (symmetric && !A->structurally_symmetric) { 61 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 62 } else if (oshift == 1) { 63 PetscInt nz = a->i[A->rmap->n]; 64 /* malloc space and add 1 to i and j indices */ 65 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 66 for (i=0; i<A->rmap->n+1; i++) (*ia)[i] = a->i[i] + 1; 67 if (ja) { 68 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 69 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 70 } 71 } else { 72 *ia = a->i; 73 if (ja) *ja = a->j; 74 } 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 80 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 81 { 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 if (!ia) PetscFunctionReturn(0); 86 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 87 ierr = PetscFree(*ia);CHKERRQ(ierr); 88 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 89 } 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 95 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 96 { 97 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 98 PetscErrorCode ierr; 99 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 100 PetscInt nz = a->i[m],row,*jj,mr,col; 101 102 PetscFunctionBegin; 103 *nn = n; 104 if (!ia) PetscFunctionReturn(0); 105 if (symmetric) { 106 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 107 } else { 108 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 109 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 110 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 111 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 112 jj = a->j; 113 for (i=0; i<nz; i++) { 114 collengths[jj[i]]++; 115 } 116 cia[0] = oshift; 117 for (i=0; i<n; i++) { 118 cia[i+1] = cia[i] + collengths[i]; 119 } 120 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 121 jj = a->j; 122 for (row=0; row<m; row++) { 123 mr = a->i[row+1] - a->i[row]; 124 for (i=0; i<mr; i++) { 125 col = *jj++; 126 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 127 } 128 } 129 ierr = PetscFree(collengths);CHKERRQ(ierr); 130 *ia = cia; *ja = cja; 131 } 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 137 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 138 { 139 PetscErrorCode ierr; 140 141 PetscFunctionBegin; 142 if (!ia) PetscFunctionReturn(0); 143 144 ierr = PetscFree(*ia);CHKERRQ(ierr); 145 ierr = PetscFree(*ja);CHKERRQ(ierr); 146 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 152 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 153 { 154 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 155 PetscInt *ai = a->i; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatSetValues_SeqAIJ" 165 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 166 { 167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 168 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 169 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 170 PetscErrorCode ierr; 171 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 172 MatScalar *ap,value,*aa = a->a; 173 PetscTruth ignorezeroentries = a->ignorezeroentries; 174 PetscTruth roworiented = a->roworiented; 175 176 PetscFunctionBegin; 177 if (v) PetscValidScalarPointer(v,6); 178 for (k=0; k<m; k++) { /* loop over added rows */ 179 row = im[k]; 180 if (row < 0) continue; 181 #if defined(PETSC_USE_DEBUG) 182 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 183 #endif 184 rp = aj + ai[row]; ap = aa + ai[row]; 185 rmax = imax[row]; nrow = ailen[row]; 186 low = 0; 187 high = nrow; 188 for (l=0; l<n; l++) { /* loop over added columns */ 189 if (in[l] < 0) continue; 190 #if defined(PETSC_USE_DEBUG) 191 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 192 #endif 193 col = in[l]; 194 if (v) { 195 if (roworiented) { 196 value = v[l + k*n]; 197 } else { 198 value = v[k + l*m]; 199 } 200 } else { 201 value = 0.; 202 } 203 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 204 205 if (col <= lastcol) low = 0; else high = nrow; 206 lastcol = col; 207 while (high-low > 5) { 208 t = (low+high)/2; 209 if (rp[t] > col) high = t; 210 else low = t; 211 } 212 for (i=low; i<high; i++) { 213 if (rp[i] > col) break; 214 if (rp[i] == col) { 215 if (is == ADD_VALUES) ap[i] += value; 216 else ap[i] = value; 217 low = i + 1; 218 goto noinsert; 219 } 220 } 221 if (value == 0.0 && ignorezeroentries) goto noinsert; 222 if (nonew == 1) goto noinsert; 223 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 224 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 225 N = nrow++ - 1; a->nz++; high++; 226 /* shift up all the later entries in this row */ 227 for (ii=N; ii>=i; ii--) { 228 rp[ii+1] = rp[ii]; 229 ap[ii+1] = ap[ii]; 230 } 231 rp[i] = col; 232 ap[i] = value; 233 low = i + 1; 234 noinsert:; 235 } 236 ailen[row] = nrow; 237 } 238 A->same_nonzero = PETSC_FALSE; 239 PetscFunctionReturn(0); 240 } 241 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "MatGetValues_SeqAIJ" 245 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 246 { 247 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 248 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 249 PetscInt *ai = a->i,*ailen = a->ilen; 250 MatScalar *ap,*aa = a->a; 251 252 PetscFunctionBegin; 253 for (k=0; k<m; k++) { /* loop over rows */ 254 row = im[k]; 255 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 256 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 257 rp = aj + ai[row]; ap = aa + ai[row]; 258 nrow = ailen[row]; 259 for (l=0; l<n; l++) { /* loop over columns */ 260 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 261 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 262 col = in[l] ; 263 high = nrow; low = 0; /* assume unsorted */ 264 while (high-low > 5) { 265 t = (low+high)/2; 266 if (rp[t] > col) high = t; 267 else low = t; 268 } 269 for (i=low; i<high; i++) { 270 if (rp[i] > col) break; 271 if (rp[i] == col) { 272 *v++ = ap[i]; 273 goto finished; 274 } 275 } 276 *v++ = 0.0; 277 finished:; 278 } 279 } 280 PetscFunctionReturn(0); 281 } 282 283 284 #undef __FUNCT__ 285 #define __FUNCT__ "MatView_SeqAIJ_Binary" 286 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 287 { 288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 289 PetscErrorCode ierr; 290 PetscInt i,*col_lens; 291 int fd; 292 293 PetscFunctionBegin; 294 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 295 ierr = PetscMalloc((4+A->rmap->n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 296 col_lens[0] = MAT_FILE_CLASSID; 297 col_lens[1] = A->rmap->n; 298 col_lens[2] = A->cmap->n; 299 col_lens[3] = a->nz; 300 301 /* store lengths of each row and write (including header) to file */ 302 for (i=0; i<A->rmap->n; i++) { 303 col_lens[4+i] = a->i[i+1] - a->i[i]; 304 } 305 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 306 ierr = PetscFree(col_lens);CHKERRQ(ierr); 307 308 /* store column indices (zero start index) */ 309 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 310 311 /* store nonzero values */ 312 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 320 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 321 { 322 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 323 PetscErrorCode ierr; 324 PetscInt i,j,m = A->rmap->n,shift=0; 325 const char *name; 326 PetscViewerFormat format; 327 328 PetscFunctionBegin; 329 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 330 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 331 if (format == PETSC_VIEWER_ASCII_MATLAB) { 332 PetscInt nofinalvalue = 0; 333 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-!shift)) { 334 nofinalvalue = 1; 335 } 336 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 337 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 338 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 339 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 340 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 341 342 for (i=0; i<m; i++) { 343 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 344 #if defined(PETSC_USE_COMPLEX) 345 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); 346 #else 347 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 348 #endif 349 } 350 } 351 if (nofinalvalue) { 352 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 353 } 354 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 355 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 356 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 357 PetscFunctionReturn(0); 358 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 359 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 360 for (i=0; i<m; i++) { 361 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 362 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 363 #if defined(PETSC_USE_COMPLEX) 364 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 365 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 366 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 367 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 368 } else if (PetscRealPart(a->a[j]) != 0.0) { 369 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 370 } 371 #else 372 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 373 #endif 374 } 375 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 376 } 377 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 378 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 379 PetscInt nzd=0,fshift=1,*sptr; 380 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 381 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 382 for (i=0; i<m; i++) { 383 sptr[i] = nzd+1; 384 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 385 if (a->j[j] >= i) { 386 #if defined(PETSC_USE_COMPLEX) 387 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 388 #else 389 if (a->a[j] != 0.0) nzd++; 390 #endif 391 } 392 } 393 } 394 sptr[m] = nzd+1; 395 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 396 for (i=0; i<m+1; i+=6) { 397 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);} 398 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);} 399 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);} 400 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 401 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 402 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 403 } 404 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 405 ierr = PetscFree(sptr);CHKERRQ(ierr); 406 for (i=0; i<m; i++) { 407 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 408 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 409 } 410 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 411 } 412 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 413 for (i=0; i<m; i++) { 414 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 415 if (a->j[j] >= i) { 416 #if defined(PETSC_USE_COMPLEX) 417 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 418 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 419 } 420 #else 421 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 422 #endif 423 } 424 } 425 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 426 } 427 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 428 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 429 PetscInt cnt = 0,jcnt; 430 PetscScalar value; 431 432 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 433 for (i=0; i<m; i++) { 434 jcnt = 0; 435 for (j=0; j<A->cmap->n; j++) { 436 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 437 value = a->a[cnt++]; 438 jcnt++; 439 } else { 440 value = 0.0; 441 } 442 #if defined(PETSC_USE_COMPLEX) 443 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 444 #else 445 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 446 #endif 447 } 448 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 449 } 450 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 451 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 452 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 453 #if defined(PETSC_USE_COMPLEX) 454 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");CHKERRQ(ierr); 455 #else 456 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");CHKERRQ(ierr); 457 #endif 458 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 459 for (i=0; i<m; i++) { 460 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 461 #if defined(PETSC_USE_COMPLEX) 462 if (PetscImaginaryPart(a->a[j]) > 0.0) { 463 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 464 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 465 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G -%G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 466 } else { 467 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 468 } 469 #else 470 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %G\n", i+shift, a->j[j]+shift, a->a[j]);CHKERRQ(ierr); 471 #endif 472 } 473 } 474 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 475 } else { 476 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 477 if (A->factortype){ 478 for (i=0; i<m; i++) { 479 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 480 /* L part */ 481 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 482 #if defined(PETSC_USE_COMPLEX) 483 if (PetscImaginaryPart(a->a[j]) > 0.0) { 484 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 485 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 486 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 487 } else { 488 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 489 } 490 #else 491 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 492 #endif 493 } 494 /* diagonal */ 495 j = a->diag[i]; 496 #if defined(PETSC_USE_COMPLEX) 497 if (PetscImaginaryPart(a->a[j]) > 0.0) { 498 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 499 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 500 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 501 } else { 502 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 503 } 504 #else 505 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 506 #endif 507 508 /* U part */ 509 for (j=a->diag[i+1]+1+shift; j<a->diag[i]+shift; j++) { 510 #if defined(PETSC_USE_COMPLEX) 511 if (PetscImaginaryPart(a->a[j]) > 0.0) { 512 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 513 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 514 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 515 } else { 516 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 517 } 518 #else 519 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 520 #endif 521 } 522 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 523 } 524 } else { 525 for (i=0; i<m; i++) { 526 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 527 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 528 #if defined(PETSC_USE_COMPLEX) 529 if (PetscImaginaryPart(a->a[j]) > 0.0) { 530 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 531 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 532 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 533 } else { 534 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 535 } 536 #else 537 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 538 #endif 539 } 540 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 541 } 542 } 543 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 544 } 545 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 #undef __FUNCT__ 550 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 551 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 552 { 553 Mat A = (Mat) Aa; 554 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 555 PetscErrorCode ierr; 556 PetscInt i,j,m = A->rmap->n,color; 557 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 558 PetscViewer viewer; 559 PetscViewerFormat format; 560 561 PetscFunctionBegin; 562 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 563 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 564 565 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 566 /* loop over matrix elements drawing boxes */ 567 568 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 569 /* Blue for negative, Cyan for zero and Red for positive */ 570 color = PETSC_DRAW_BLUE; 571 for (i=0; i<m; i++) { 572 y_l = m - i - 1.0; y_r = y_l + 1.0; 573 for (j=a->i[i]; j<a->i[i+1]; j++) { 574 x_l = a->j[j] ; x_r = x_l + 1.0; 575 #if defined(PETSC_USE_COMPLEX) 576 if (PetscRealPart(a->a[j]) >= 0.) continue; 577 #else 578 if (a->a[j] >= 0.) continue; 579 #endif 580 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 581 } 582 } 583 color = PETSC_DRAW_CYAN; 584 for (i=0; i<m; i++) { 585 y_l = m - i - 1.0; y_r = y_l + 1.0; 586 for (j=a->i[i]; j<a->i[i+1]; j++) { 587 x_l = a->j[j]; x_r = x_l + 1.0; 588 if (a->a[j] != 0.) continue; 589 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 590 } 591 } 592 color = PETSC_DRAW_RED; 593 for (i=0; i<m; i++) { 594 y_l = m - i - 1.0; y_r = y_l + 1.0; 595 for (j=a->i[i]; j<a->i[i+1]; j++) { 596 x_l = a->j[j]; x_r = x_l + 1.0; 597 #if defined(PETSC_USE_COMPLEX) 598 if (PetscRealPart(a->a[j]) <= 0.) continue; 599 #else 600 if (a->a[j] <= 0.) continue; 601 #endif 602 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 603 } 604 } 605 } else { 606 /* use contour shading to indicate magnitude of values */ 607 /* first determine max of all nonzero values */ 608 PetscInt nz = a->nz,count; 609 PetscDraw popup; 610 PetscReal scale; 611 612 for (i=0; i<nz; i++) { 613 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 614 } 615 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 616 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 617 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 618 count = 0; 619 for (i=0; i<m; i++) { 620 y_l = m - i - 1.0; y_r = y_l + 1.0; 621 for (j=a->i[i]; j<a->i[i+1]; j++) { 622 x_l = a->j[j]; x_r = x_l + 1.0; 623 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 624 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 625 count++; 626 } 627 } 628 } 629 PetscFunctionReturn(0); 630 } 631 632 #undef __FUNCT__ 633 #define __FUNCT__ "MatView_SeqAIJ_Draw" 634 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 635 { 636 PetscErrorCode ierr; 637 PetscDraw draw; 638 PetscReal xr,yr,xl,yl,h,w; 639 PetscTruth isnull; 640 641 PetscFunctionBegin; 642 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 643 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 644 if (isnull) PetscFunctionReturn(0); 645 646 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 647 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 648 xr += w; yr += h; xl = -w; yl = -h; 649 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 650 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 651 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 652 PetscFunctionReturn(0); 653 } 654 655 #undef __FUNCT__ 656 #define __FUNCT__ "MatView_SeqAIJ" 657 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 658 { 659 PetscErrorCode ierr; 660 PetscTruth iascii,isbinary,isdraw; 661 662 PetscFunctionBegin; 663 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 664 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 665 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 666 if (iascii) { 667 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 668 } else if (isbinary) { 669 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 670 } else if (isdraw) { 671 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 672 } else { 673 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 674 } 675 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 679 #undef __FUNCT__ 680 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 681 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 682 { 683 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 684 PetscErrorCode ierr; 685 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 686 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 687 MatScalar *aa = a->a,*ap; 688 PetscReal ratio=0.6; 689 690 PetscFunctionBegin; 691 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 692 693 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 694 for (i=1; i<m; i++) { 695 /* move each row back by the amount of empty slots (fshift) before it*/ 696 fshift += imax[i-1] - ailen[i-1]; 697 rmax = PetscMax(rmax,ailen[i]); 698 if (fshift) { 699 ip = aj + ai[i] ; 700 ap = aa + ai[i] ; 701 N = ailen[i]; 702 for (j=0; j<N; j++) { 703 ip[j-fshift] = ip[j]; 704 ap[j-fshift] = ap[j]; 705 } 706 } 707 ai[i] = ai[i-1] + ailen[i-1]; 708 } 709 if (m) { 710 fshift += imax[m-1] - ailen[m-1]; 711 ai[m] = ai[m-1] + ailen[m-1]; 712 } 713 /* reset ilen and imax for each row */ 714 for (i=0; i<m; i++) { 715 ailen[i] = imax[i] = ai[i+1] - ai[i]; 716 } 717 a->nz = ai[m]; 718 if (fshift && a->nounused == -1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D, %D unneeded", m, A->cmap->n, fshift); 719 720 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 721 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 722 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 723 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 724 A->info.mallocs += a->reallocs; 725 a->reallocs = 0; 726 A->info.nz_unneeded = (double)fshift; 727 a->rmax = rmax; 728 729 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 730 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 731 A->same_nonzero = PETSC_TRUE; 732 733 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 734 735 a->idiagvalid = PETSC_FALSE; 736 PetscFunctionReturn(0); 737 } 738 739 #undef __FUNCT__ 740 #define __FUNCT__ "MatRealPart_SeqAIJ" 741 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 742 { 743 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 744 PetscInt i,nz = a->nz; 745 MatScalar *aa = a->a; 746 747 PetscFunctionBegin; 748 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 749 PetscFunctionReturn(0); 750 } 751 752 #undef __FUNCT__ 753 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 754 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 755 { 756 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 757 PetscInt i,nz = a->nz; 758 MatScalar *aa = a->a; 759 760 PetscFunctionBegin; 761 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 762 PetscFunctionReturn(0); 763 } 764 765 #undef __FUNCT__ 766 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 767 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 768 { 769 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 770 PetscErrorCode ierr; 771 772 PetscFunctionBegin; 773 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 774 PetscFunctionReturn(0); 775 } 776 777 #undef __FUNCT__ 778 #define __FUNCT__ "MatDestroy_SeqAIJ" 779 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 780 { 781 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 782 PetscErrorCode ierr; 783 784 PetscFunctionBegin; 785 #if defined(PETSC_USE_LOG) 786 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 787 #endif 788 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 789 if (a->row) { 790 ierr = ISDestroy(a->row);CHKERRQ(ierr); 791 } 792 if (a->col) { 793 ierr = ISDestroy(a->col);CHKERRQ(ierr); 794 } 795 ierr = PetscFree(a->diag);CHKERRQ(ierr); 796 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 797 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 798 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 799 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 800 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 801 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 802 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 803 if (a->XtoY) {ierr = MatDestroy(a->XtoY);CHKERRQ(ierr);} 804 if (a->compressedrow.checked && a->compressedrow.use){ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);} 805 806 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 807 808 ierr = PetscFree(a);CHKERRQ(ierr); 809 810 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 811 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 812 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 813 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 814 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 815 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 816 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqcsrperm_C","",PETSC_NULL);CHKERRQ(ierr); 817 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 818 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 819 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 820 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 821 PetscFunctionReturn(0); 822 } 823 824 #undef __FUNCT__ 825 #define __FUNCT__ "MatSetOption_SeqAIJ" 826 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscTruth flg) 827 { 828 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 829 PetscErrorCode ierr; 830 831 PetscFunctionBegin; 832 switch (op) { 833 case MAT_ROW_ORIENTED: 834 a->roworiented = flg; 835 break; 836 case MAT_KEEP_NONZERO_PATTERN: 837 a->keepnonzeropattern = flg; 838 break; 839 case MAT_NEW_NONZERO_LOCATIONS: 840 a->nonew = (flg ? 0 : 1); 841 break; 842 case MAT_NEW_NONZERO_LOCATION_ERR: 843 a->nonew = (flg ? -1 : 0); 844 break; 845 case MAT_NEW_NONZERO_ALLOCATION_ERR: 846 a->nonew = (flg ? -2 : 0); 847 break; 848 case MAT_UNUSED_NONZERO_LOCATION_ERR: 849 a->nounused = (flg ? -1 : 0); 850 break; 851 case MAT_IGNORE_ZERO_ENTRIES: 852 a->ignorezeroentries = flg; 853 break; 854 case MAT_USE_COMPRESSEDROW: 855 a->compressedrow.use = flg; 856 break; 857 case MAT_SYMMETRIC: 858 case MAT_STRUCTURALLY_SYMMETRIC: 859 case MAT_HERMITIAN: 860 case MAT_SYMMETRY_ETERNAL: 861 case MAT_NEW_DIAGONALS: 862 case MAT_IGNORE_OFF_PROC_ENTRIES: 863 case MAT_USE_HASH_TABLE: 864 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 865 break; 866 case MAT_USE_INODES: 867 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 868 break; 869 default: 870 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 871 } 872 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 873 PetscFunctionReturn(0); 874 } 875 876 #undef __FUNCT__ 877 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 878 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 879 { 880 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 881 PetscErrorCode ierr; 882 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 883 PetscScalar *aa=a->a,*x,zero=0.0; 884 885 PetscFunctionBegin; 886 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 887 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 888 889 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU){ 890 PetscInt *diag=a->diag; 891 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 892 for (i=0; i<n; i++) x[i] = aa[diag[i]]; 893 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 894 PetscFunctionReturn(0); 895 } 896 897 ierr = VecSet(v,zero);CHKERRQ(ierr); 898 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 899 for (i=0; i<n; i++) { 900 nz = ai[i+1] - ai[i]; 901 if (!nz) x[i] = 0.0; 902 for (j=ai[i]; j<ai[i+1]; j++){ 903 if (aj[j] == i) { 904 x[i] = aa[j]; 905 break; 906 } 907 } 908 } 909 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 910 PetscFunctionReturn(0); 911 } 912 913 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 914 #undef __FUNCT__ 915 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 916 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 917 { 918 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 919 PetscScalar *x,*y; 920 PetscErrorCode ierr; 921 PetscInt m = A->rmap->n; 922 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 923 MatScalar *v; 924 PetscScalar alpha; 925 PetscInt n,i,j,*idx,*ii,*ridx=PETSC_NULL; 926 Mat_CompressedRow cprow = a->compressedrow; 927 PetscTruth usecprow = cprow.use; 928 #endif 929 930 PetscFunctionBegin; 931 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 932 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 933 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 934 935 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 936 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 937 #else 938 if (usecprow){ 939 m = cprow.nrows; 940 ii = cprow.i; 941 ridx = cprow.rindex; 942 } else { 943 ii = a->i; 944 } 945 for (i=0; i<m; i++) { 946 idx = a->j + ii[i] ; 947 v = a->a + ii[i] ; 948 n = ii[i+1] - ii[i]; 949 if (usecprow){ 950 alpha = x[ridx[i]]; 951 } else { 952 alpha = x[i]; 953 } 954 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 955 } 956 #endif 957 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 958 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 959 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 960 PetscFunctionReturn(0); 961 } 962 963 #undef __FUNCT__ 964 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 965 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 966 { 967 PetscErrorCode ierr; 968 969 PetscFunctionBegin; 970 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 971 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 972 PetscFunctionReturn(0); 973 } 974 975 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 976 #undef __FUNCT__ 977 #define __FUNCT__ "MatMult_SeqAIJ" 978 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 979 { 980 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 981 PetscScalar *y; 982 const PetscScalar *x; 983 const MatScalar *aa; 984 PetscErrorCode ierr; 985 PetscInt m=A->rmap->n; 986 const PetscInt *aj,*ii,*ridx=PETSC_NULL; 987 PetscInt n,i,nonzerorow=0; 988 PetscScalar sum; 989 PetscTruth usecprow=a->compressedrow.use; 990 991 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 992 #pragma disjoint(*x,*y,*aa) 993 #endif 994 995 PetscFunctionBegin; 996 ierr = VecGetArray(xx,(PetscScalar**)&x);CHKERRQ(ierr); 997 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 998 aj = a->j; 999 aa = a->a; 1000 ii = a->i; 1001 if (usecprow){ /* use compressed row format */ 1002 m = a->compressedrow.nrows; 1003 ii = a->compressedrow.i; 1004 ridx = a->compressedrow.rindex; 1005 for (i=0; i<m; i++){ 1006 n = ii[i+1] - ii[i]; 1007 aj = a->j + ii[i]; 1008 aa = a->a + ii[i]; 1009 sum = 0.0; 1010 nonzerorow += (n>0); 1011 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1012 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1013 y[*ridx++] = sum; 1014 } 1015 } else { /* do not use compressed row format */ 1016 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1017 fortranmultaij_(&m,x,ii,aj,aa,y); 1018 #else 1019 for (i=0; i<m; i++) { 1020 n = ii[i+1] - ii[i]; 1021 aj = a->j + ii[i]; 1022 aa = a->a + ii[i]; 1023 sum = 0.0; 1024 nonzerorow += (n>0); 1025 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1026 y[i] = sum; 1027 } 1028 #endif 1029 } 1030 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1031 ierr = VecRestoreArray(xx,(PetscScalar**)&x);CHKERRQ(ierr); 1032 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1033 PetscFunctionReturn(0); 1034 } 1035 1036 #include "../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h" 1037 #undef __FUNCT__ 1038 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1039 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1040 { 1041 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1042 PetscScalar *x,*y,*z; 1043 const MatScalar *aa; 1044 PetscErrorCode ierr; 1045 PetscInt m = A->rmap->n,*aj,*ii; 1046 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1047 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 1048 PetscScalar sum; 1049 PetscTruth usecprow=a->compressedrow.use; 1050 #endif 1051 1052 PetscFunctionBegin; 1053 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1054 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1055 if (zz != yy) { 1056 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1057 } else { 1058 z = y; 1059 } 1060 1061 aj = a->j; 1062 aa = a->a; 1063 ii = a->i; 1064 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1065 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1066 #else 1067 if (usecprow){ /* use compressed row format */ 1068 if (zz != yy){ 1069 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1070 } 1071 m = a->compressedrow.nrows; 1072 ii = a->compressedrow.i; 1073 ridx = a->compressedrow.rindex; 1074 for (i=0; i<m; i++){ 1075 n = ii[i+1] - ii[i]; 1076 aj = a->j + ii[i]; 1077 aa = a->a + ii[i]; 1078 sum = y[*ridx]; 1079 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1080 z[*ridx++] = sum; 1081 } 1082 } else { /* do not use compressed row format */ 1083 for (i=0; i<m; i++) { 1084 jrow = ii[i]; 1085 n = ii[i+1] - jrow; 1086 sum = y[i]; 1087 for (j=0; j<n; j++) { 1088 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1089 } 1090 z[i] = sum; 1091 } 1092 } 1093 #endif 1094 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1095 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1096 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1097 if (zz != yy) { 1098 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1099 } 1100 PetscFunctionReturn(0); 1101 } 1102 1103 /* 1104 Adds diagonal pointers to sparse matrix structure. 1105 */ 1106 #undef __FUNCT__ 1107 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1108 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1109 { 1110 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1111 PetscErrorCode ierr; 1112 PetscInt i,j,m = A->rmap->n; 1113 1114 PetscFunctionBegin; 1115 if (!a->diag) { 1116 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1117 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1118 } 1119 for (i=0; i<A->rmap->n; i++) { 1120 a->diag[i] = a->i[i+1]; 1121 for (j=a->i[i]; j<a->i[i+1]; j++) { 1122 if (a->j[j] == i) { 1123 a->diag[i] = j; 1124 break; 1125 } 1126 } 1127 } 1128 PetscFunctionReturn(0); 1129 } 1130 1131 /* 1132 Checks for missing diagonals 1133 */ 1134 #undef __FUNCT__ 1135 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1136 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1137 { 1138 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1139 PetscInt *diag,*jj = a->j,i; 1140 1141 PetscFunctionBegin; 1142 *missing = PETSC_FALSE; 1143 if (A->rmap->n > 0 && !jj) { 1144 *missing = PETSC_TRUE; 1145 if (d) *d = 0; 1146 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1147 } else { 1148 diag = a->diag; 1149 for (i=0; i<A->rmap->n; i++) { 1150 if (jj[diag[i]] != i) { 1151 *missing = PETSC_TRUE; 1152 if (d) *d = i; 1153 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1154 } 1155 } 1156 } 1157 PetscFunctionReturn(0); 1158 } 1159 1160 EXTERN_C_BEGIN 1161 #undef __FUNCT__ 1162 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1163 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1164 { 1165 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1166 PetscErrorCode ierr; 1167 PetscInt i,*diag,m = A->rmap->n; 1168 MatScalar *v = a->a; 1169 PetscScalar *idiag,*mdiag; 1170 1171 PetscFunctionBegin; 1172 if (a->idiagvalid) PetscFunctionReturn(0); 1173 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1174 diag = a->diag; 1175 if (!a->idiag) { 1176 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1177 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1178 v = a->a; 1179 } 1180 mdiag = a->mdiag; 1181 idiag = a->idiag; 1182 1183 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1184 for (i=0; i<m; i++) { 1185 mdiag[i] = v[diag[i]]; 1186 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1187 idiag[i] = 1.0/v[diag[i]]; 1188 } 1189 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1190 } else { 1191 for (i=0; i<m; i++) { 1192 mdiag[i] = v[diag[i]]; 1193 idiag[i] = omega/(fshift + v[diag[i]]); 1194 } 1195 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1196 } 1197 a->idiagvalid = PETSC_TRUE; 1198 PetscFunctionReturn(0); 1199 } 1200 EXTERN_C_END 1201 1202 #include "../src/mat/impls/aij/seq/ftn-kernels/frelax.h" 1203 #undef __FUNCT__ 1204 #define __FUNCT__ "MatSOR_SeqAIJ" 1205 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1206 { 1207 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1208 PetscScalar *x,d,sum,*t,scale; 1209 const MatScalar *v = a->a,*idiag=0,*mdiag; 1210 const PetscScalar *b, *bs,*xb, *ts; 1211 PetscErrorCode ierr; 1212 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1213 const PetscInt *idx,*diag; 1214 1215 PetscFunctionBegin; 1216 its = its*lits; 1217 1218 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1219 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1220 a->fshift = fshift; 1221 a->omega = omega; 1222 1223 diag = a->diag; 1224 t = a->ssor_work; 1225 idiag = a->idiag; 1226 mdiag = a->mdiag; 1227 1228 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1229 if (xx != bb) { 1230 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1231 } else { 1232 b = x; 1233 } 1234 CHKMEMQ; 1235 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1236 if (flag == SOR_APPLY_UPPER) { 1237 /* apply (U + D/omega) to the vector */ 1238 bs = b; 1239 for (i=0; i<m; i++) { 1240 d = fshift + mdiag[i]; 1241 n = a->i[i+1] - diag[i] - 1; 1242 idx = a->j + diag[i] + 1; 1243 v = a->a + diag[i] + 1; 1244 sum = b[i]*d/omega; 1245 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1246 x[i] = sum; 1247 } 1248 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1249 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1250 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1251 PetscFunctionReturn(0); 1252 } 1253 1254 if (flag == SOR_APPLY_LOWER) { 1255 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1256 } else if (flag & SOR_EISENSTAT) { 1257 /* Let A = L + U + D; where L is lower trianglar, 1258 U is upper triangular, E = D/omega; This routine applies 1259 1260 (L + E)^{-1} A (U + E)^{-1} 1261 1262 to a vector efficiently using Eisenstat's trick. 1263 */ 1264 scale = (2.0/omega) - 1.0; 1265 1266 /* x = (E + U)^{-1} b */ 1267 for (i=m-1; i>=0; i--) { 1268 n = a->i[i+1] - diag[i] - 1; 1269 idx = a->j + diag[i] + 1; 1270 v = a->a + diag[i] + 1; 1271 sum = b[i]; 1272 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1273 x[i] = sum*idiag[i]; 1274 } 1275 1276 /* t = b - (2*E - D)x */ 1277 v = a->a; 1278 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1279 1280 /* t = (E + L)^{-1}t */ 1281 ts = t; 1282 diag = a->diag; 1283 for (i=0; i<m; i++) { 1284 n = diag[i] - a->i[i]; 1285 idx = a->j + a->i[i]; 1286 v = a->a + a->i[i]; 1287 sum = t[i]; 1288 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1289 t[i] = sum*idiag[i]; 1290 /* x = x + t */ 1291 x[i] += t[i]; 1292 } 1293 1294 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1295 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1296 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1297 PetscFunctionReturn(0); 1298 } 1299 if (flag & SOR_ZERO_INITIAL_GUESS) { 1300 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1301 for (i=0; i<m; i++) { 1302 n = diag[i] - a->i[i]; 1303 idx = a->j + a->i[i]; 1304 v = a->a + a->i[i]; 1305 sum = b[i]; 1306 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1307 t[i] = sum; 1308 x[i] = sum*idiag[i]; 1309 } 1310 xb = t; 1311 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1312 } else xb = b; 1313 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1314 for (i=m-1; i>=0; i--) { 1315 n = a->i[i+1] - diag[i] - 1; 1316 idx = a->j + diag[i] + 1; 1317 v = a->a + diag[i] + 1; 1318 sum = xb[i]; 1319 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1320 if (xb == b) { 1321 x[i] = sum*idiag[i]; 1322 } else { 1323 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1324 } 1325 } 1326 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1327 } 1328 its--; 1329 } 1330 while (its--) { 1331 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1332 for (i=0; i<m; i++) { 1333 n = a->i[i+1] - a->i[i]; 1334 idx = a->j + a->i[i]; 1335 v = a->a + a->i[i]; 1336 sum = b[i]; 1337 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1338 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1339 } 1340 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1341 } 1342 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1343 for (i=m-1; i>=0; i--) { 1344 n = a->i[i+1] - a->i[i]; 1345 idx = a->j + a->i[i]; 1346 v = a->a + a->i[i]; 1347 sum = b[i]; 1348 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1349 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1350 } 1351 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1352 } 1353 } 1354 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1355 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1356 CHKMEMQ; PetscFunctionReturn(0); 1357 } 1358 1359 1360 #undef __FUNCT__ 1361 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1362 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1363 { 1364 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1365 1366 PetscFunctionBegin; 1367 info->block_size = 1.0; 1368 info->nz_allocated = (double)a->maxnz; 1369 info->nz_used = (double)a->nz; 1370 info->nz_unneeded = (double)(a->maxnz - a->nz); 1371 info->assemblies = (double)A->num_ass; 1372 info->mallocs = (double)A->info.mallocs; 1373 info->memory = ((PetscObject)A)->mem; 1374 if (A->factortype) { 1375 info->fill_ratio_given = A->info.fill_ratio_given; 1376 info->fill_ratio_needed = A->info.fill_ratio_needed; 1377 info->factor_mallocs = A->info.factor_mallocs; 1378 } else { 1379 info->fill_ratio_given = 0; 1380 info->fill_ratio_needed = 0; 1381 info->factor_mallocs = 0; 1382 } 1383 PetscFunctionReturn(0); 1384 } 1385 1386 #undef __FUNCT__ 1387 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1388 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1389 { 1390 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1391 PetscInt i,m = A->rmap->n - 1,d = 0; 1392 PetscErrorCode ierr; 1393 PetscTruth missing; 1394 1395 PetscFunctionBegin; 1396 if (a->keepnonzeropattern) { 1397 for (i=0; i<N; i++) { 1398 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1399 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1400 } 1401 if (diag != 0.0) { 1402 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1403 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1404 for (i=0; i<N; i++) { 1405 a->a[a->diag[rows[i]]] = diag; 1406 } 1407 } 1408 A->same_nonzero = PETSC_TRUE; 1409 } else { 1410 if (diag != 0.0) { 1411 for (i=0; i<N; i++) { 1412 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1413 if (a->ilen[rows[i]] > 0) { 1414 a->ilen[rows[i]] = 1; 1415 a->a[a->i[rows[i]]] = diag; 1416 a->j[a->i[rows[i]]] = rows[i]; 1417 } else { /* in case row was completely empty */ 1418 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1419 } 1420 } 1421 } else { 1422 for (i=0; i<N; i++) { 1423 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1424 a->ilen[rows[i]] = 0; 1425 } 1426 } 1427 A->same_nonzero = PETSC_FALSE; 1428 } 1429 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1430 PetscFunctionReturn(0); 1431 } 1432 1433 #undef __FUNCT__ 1434 #define __FUNCT__ "MatGetRow_SeqAIJ" 1435 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1436 { 1437 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1438 PetscInt *itmp; 1439 1440 PetscFunctionBegin; 1441 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1442 1443 *nz = a->i[row+1] - a->i[row]; 1444 if (v) *v = a->a + a->i[row]; 1445 if (idx) { 1446 itmp = a->j + a->i[row]; 1447 if (*nz) { 1448 *idx = itmp; 1449 } 1450 else *idx = 0; 1451 } 1452 PetscFunctionReturn(0); 1453 } 1454 1455 /* remove this function? */ 1456 #undef __FUNCT__ 1457 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1458 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1459 { 1460 PetscFunctionBegin; 1461 PetscFunctionReturn(0); 1462 } 1463 1464 #undef __FUNCT__ 1465 #define __FUNCT__ "MatNorm_SeqAIJ" 1466 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1467 { 1468 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1469 MatScalar *v = a->a; 1470 PetscReal sum = 0.0; 1471 PetscErrorCode ierr; 1472 PetscInt i,j; 1473 1474 PetscFunctionBegin; 1475 if (type == NORM_FROBENIUS) { 1476 for (i=0; i<a->nz; i++) { 1477 #if defined(PETSC_USE_COMPLEX) 1478 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1479 #else 1480 sum += (*v)*(*v); v++; 1481 #endif 1482 } 1483 *nrm = sqrt(sum); 1484 } else if (type == NORM_1) { 1485 PetscReal *tmp; 1486 PetscInt *jj = a->j; 1487 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1488 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1489 *nrm = 0.0; 1490 for (j=0; j<a->nz; j++) { 1491 tmp[*jj++] += PetscAbsScalar(*v); v++; 1492 } 1493 for (j=0; j<A->cmap->n; j++) { 1494 if (tmp[j] > *nrm) *nrm = tmp[j]; 1495 } 1496 ierr = PetscFree(tmp);CHKERRQ(ierr); 1497 } else if (type == NORM_INFINITY) { 1498 *nrm = 0.0; 1499 for (j=0; j<A->rmap->n; j++) { 1500 v = a->a + a->i[j]; 1501 sum = 0.0; 1502 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1503 sum += PetscAbsScalar(*v); v++; 1504 } 1505 if (sum > *nrm) *nrm = sum; 1506 } 1507 } else { 1508 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1509 } 1510 PetscFunctionReturn(0); 1511 } 1512 1513 #undef __FUNCT__ 1514 #define __FUNCT__ "MatTranspose_SeqAIJ" 1515 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1516 { 1517 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1518 Mat C; 1519 PetscErrorCode ierr; 1520 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1521 MatScalar *array = a->a; 1522 1523 PetscFunctionBegin; 1524 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1525 1526 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1527 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1528 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1529 1530 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1531 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1532 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1533 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1534 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1535 ierr = PetscFree(col);CHKERRQ(ierr); 1536 } else { 1537 C = *B; 1538 } 1539 1540 for (i=0; i<m; i++) { 1541 len = ai[i+1]-ai[i]; 1542 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1543 array += len; 1544 aj += len; 1545 } 1546 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1547 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1548 1549 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1550 *B = C; 1551 } else { 1552 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1553 } 1554 PetscFunctionReturn(0); 1555 } 1556 1557 EXTERN_C_BEGIN 1558 #undef __FUNCT__ 1559 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1560 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1561 { 1562 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1563 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1564 MatScalar *va,*vb; 1565 PetscErrorCode ierr; 1566 PetscInt ma,na,mb,nb, i; 1567 1568 PetscFunctionBegin; 1569 bij = (Mat_SeqAIJ *) B->data; 1570 1571 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1572 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1573 if (ma!=nb || na!=mb){ 1574 *f = PETSC_FALSE; 1575 PetscFunctionReturn(0); 1576 } 1577 aii = aij->i; bii = bij->i; 1578 adx = aij->j; bdx = bij->j; 1579 va = aij->a; vb = bij->a; 1580 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1581 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1582 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1583 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1584 1585 *f = PETSC_TRUE; 1586 for (i=0; i<ma; i++) { 1587 while (aptr[i]<aii[i+1]) { 1588 PetscInt idc,idr; 1589 PetscScalar vc,vr; 1590 /* column/row index/value */ 1591 idc = adx[aptr[i]]; 1592 idr = bdx[bptr[idc]]; 1593 vc = va[aptr[i]]; 1594 vr = vb[bptr[idc]]; 1595 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1596 *f = PETSC_FALSE; 1597 goto done; 1598 } else { 1599 aptr[i]++; 1600 if (B || i!=idc) bptr[idc]++; 1601 } 1602 } 1603 } 1604 done: 1605 ierr = PetscFree(aptr);CHKERRQ(ierr); 1606 if (B) { 1607 ierr = PetscFree(bptr);CHKERRQ(ierr); 1608 } 1609 PetscFunctionReturn(0); 1610 } 1611 EXTERN_C_END 1612 1613 EXTERN_C_BEGIN 1614 #undef __FUNCT__ 1615 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1616 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1617 { 1618 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1619 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1620 MatScalar *va,*vb; 1621 PetscErrorCode ierr; 1622 PetscInt ma,na,mb,nb, i; 1623 1624 PetscFunctionBegin; 1625 bij = (Mat_SeqAIJ *) B->data; 1626 1627 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1628 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1629 if (ma!=nb || na!=mb){ 1630 *f = PETSC_FALSE; 1631 PetscFunctionReturn(0); 1632 } 1633 aii = aij->i; bii = bij->i; 1634 adx = aij->j; bdx = bij->j; 1635 va = aij->a; vb = bij->a; 1636 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1637 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1638 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1639 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1640 1641 *f = PETSC_TRUE; 1642 for (i=0; i<ma; i++) { 1643 while (aptr[i]<aii[i+1]) { 1644 PetscInt idc,idr; 1645 PetscScalar vc,vr; 1646 /* column/row index/value */ 1647 idc = adx[aptr[i]]; 1648 idr = bdx[bptr[idc]]; 1649 vc = va[aptr[i]]; 1650 vr = vb[bptr[idc]]; 1651 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1652 *f = PETSC_FALSE; 1653 goto done; 1654 } else { 1655 aptr[i]++; 1656 if (B || i!=idc) bptr[idc]++; 1657 } 1658 } 1659 } 1660 done: 1661 ierr = PetscFree(aptr);CHKERRQ(ierr); 1662 if (B) { 1663 ierr = PetscFree(bptr);CHKERRQ(ierr); 1664 } 1665 PetscFunctionReturn(0); 1666 } 1667 EXTERN_C_END 1668 1669 #undef __FUNCT__ 1670 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1671 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1672 { 1673 PetscErrorCode ierr; 1674 PetscFunctionBegin; 1675 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1676 PetscFunctionReturn(0); 1677 } 1678 1679 #undef __FUNCT__ 1680 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1681 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1682 { 1683 PetscErrorCode ierr; 1684 PetscFunctionBegin; 1685 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1686 PetscFunctionReturn(0); 1687 } 1688 1689 #undef __FUNCT__ 1690 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1691 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1692 { 1693 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1694 PetscScalar *l,*r,x; 1695 MatScalar *v; 1696 PetscErrorCode ierr; 1697 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1698 1699 PetscFunctionBegin; 1700 if (ll) { 1701 /* The local size is used so that VecMPI can be passed to this routine 1702 by MatDiagonalScale_MPIAIJ */ 1703 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1704 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1705 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1706 v = a->a; 1707 for (i=0; i<m; i++) { 1708 x = l[i]; 1709 M = a->i[i+1] - a->i[i]; 1710 for (j=0; j<M; j++) { (*v++) *= x;} 1711 } 1712 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1713 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1714 } 1715 if (rr) { 1716 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1717 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1718 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1719 v = a->a; jj = a->j; 1720 for (i=0; i<nz; i++) { 1721 (*v++) *= r[*jj++]; 1722 } 1723 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1724 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1725 } 1726 PetscFunctionReturn(0); 1727 } 1728 1729 #undef __FUNCT__ 1730 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1731 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1732 { 1733 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1734 PetscErrorCode ierr; 1735 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1736 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1737 const PetscInt *irow,*icol; 1738 PetscInt nrows,ncols; 1739 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1740 MatScalar *a_new,*mat_a; 1741 Mat C; 1742 PetscTruth stride,sorted; 1743 1744 PetscFunctionBegin; 1745 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1746 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1747 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1748 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1749 1750 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1751 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1752 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1753 1754 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1755 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1756 if (stride && step == 1) { 1757 /* special case of contiguous rows */ 1758 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1759 /* loop over new rows determining lens and starting points */ 1760 for (i=0; i<nrows; i++) { 1761 kstart = ai[irow[i]]; 1762 kend = kstart + ailen[irow[i]]; 1763 for (k=kstart; k<kend; k++) { 1764 if (aj[k] >= first) { 1765 starts[i] = k; 1766 break; 1767 } 1768 } 1769 sum = 0; 1770 while (k < kend) { 1771 if (aj[k++] >= first+ncols) break; 1772 sum++; 1773 } 1774 lens[i] = sum; 1775 } 1776 /* create submatrix */ 1777 if (scall == MAT_REUSE_MATRIX) { 1778 PetscInt n_cols,n_rows; 1779 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1780 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1781 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1782 C = *B; 1783 } else { 1784 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1785 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1786 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1787 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1788 } 1789 c = (Mat_SeqAIJ*)C->data; 1790 1791 /* loop over rows inserting into submatrix */ 1792 a_new = c->a; 1793 j_new = c->j; 1794 i_new = c->i; 1795 1796 for (i=0; i<nrows; i++) { 1797 ii = starts[i]; 1798 lensi = lens[i]; 1799 for (k=0; k<lensi; k++) { 1800 *j_new++ = aj[ii+k] - first; 1801 } 1802 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1803 a_new += lensi; 1804 i_new[i+1] = i_new[i] + lensi; 1805 c->ilen[i] = lensi; 1806 } 1807 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1808 } else { 1809 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1810 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1811 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1812 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1813 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1814 /* determine lens of each row */ 1815 for (i=0; i<nrows; i++) { 1816 kstart = ai[irow[i]]; 1817 kend = kstart + a->ilen[irow[i]]; 1818 lens[i] = 0; 1819 for (k=kstart; k<kend; k++) { 1820 if (smap[aj[k]]) { 1821 lens[i]++; 1822 } 1823 } 1824 } 1825 /* Create and fill new matrix */ 1826 if (scall == MAT_REUSE_MATRIX) { 1827 PetscTruth equal; 1828 1829 c = (Mat_SeqAIJ *)((*B)->data); 1830 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1831 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1832 if (!equal) { 1833 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1834 } 1835 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1836 C = *B; 1837 } else { 1838 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1839 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1840 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1841 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1842 } 1843 c = (Mat_SeqAIJ *)(C->data); 1844 for (i=0; i<nrows; i++) { 1845 row = irow[i]; 1846 kstart = ai[row]; 1847 kend = kstart + a->ilen[row]; 1848 mat_i = c->i[i]; 1849 mat_j = c->j + mat_i; 1850 mat_a = c->a + mat_i; 1851 mat_ilen = c->ilen + i; 1852 for (k=kstart; k<kend; k++) { 1853 if ((tcol=smap[a->j[k]])) { 1854 *mat_j++ = tcol - 1; 1855 *mat_a++ = a->a[k]; 1856 (*mat_ilen)++; 1857 1858 } 1859 } 1860 } 1861 /* Free work space */ 1862 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1863 ierr = PetscFree(smap);CHKERRQ(ierr); 1864 ierr = PetscFree(lens);CHKERRQ(ierr); 1865 } 1866 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1867 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1868 1869 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1870 *B = C; 1871 PetscFunctionReturn(0); 1872 } 1873 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1876 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1877 { 1878 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1879 PetscErrorCode ierr; 1880 Mat outA; 1881 PetscTruth row_identity,col_identity; 1882 1883 PetscFunctionBegin; 1884 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1885 1886 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1887 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1888 1889 outA = inA; 1890 outA->factortype = MAT_FACTOR_LU; 1891 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1892 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1893 a->row = row; 1894 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1895 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1896 a->col = col; 1897 1898 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1899 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1900 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1901 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1902 1903 if (!a->solve_work) { /* this matrix may have been factored before */ 1904 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1905 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1906 } 1907 1908 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1909 if (row_identity && col_identity) { 1910 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1911 } else { 1912 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1913 } 1914 PetscFunctionReturn(0); 1915 } 1916 1917 #undef __FUNCT__ 1918 #define __FUNCT__ "MatScale_SeqAIJ" 1919 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1920 { 1921 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1922 PetscScalar oalpha = alpha; 1923 PetscErrorCode ierr; 1924 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1925 1926 PetscFunctionBegin; 1927 BLASscal_(&bnz,&oalpha,a->a,&one); 1928 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1929 PetscFunctionReturn(0); 1930 } 1931 1932 #undef __FUNCT__ 1933 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1934 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1935 { 1936 PetscErrorCode ierr; 1937 PetscInt i; 1938 1939 PetscFunctionBegin; 1940 if (scall == MAT_INITIAL_MATRIX) { 1941 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1942 } 1943 1944 for (i=0; i<n; i++) { 1945 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1946 } 1947 PetscFunctionReturn(0); 1948 } 1949 1950 #undef __FUNCT__ 1951 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1952 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1953 { 1954 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1955 PetscErrorCode ierr; 1956 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1957 const PetscInt *idx; 1958 PetscInt start,end,*ai,*aj; 1959 PetscBT table; 1960 1961 PetscFunctionBegin; 1962 m = A->rmap->n; 1963 ai = a->i; 1964 aj = a->j; 1965 1966 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1967 1968 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1969 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1970 1971 for (i=0; i<is_max; i++) { 1972 /* Initialize the two local arrays */ 1973 isz = 0; 1974 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1975 1976 /* Extract the indices, assume there can be duplicate entries */ 1977 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1978 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1979 1980 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1981 for (j=0; j<n ; ++j){ 1982 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1983 } 1984 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1985 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1986 1987 k = 0; 1988 for (j=0; j<ov; j++){ /* for each overlap */ 1989 n = isz; 1990 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1991 row = nidx[k]; 1992 start = ai[row]; 1993 end = ai[row+1]; 1994 for (l = start; l<end ; l++){ 1995 val = aj[l] ; 1996 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1997 } 1998 } 1999 } 2000 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 2001 } 2002 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 2003 ierr = PetscFree(nidx);CHKERRQ(ierr); 2004 PetscFunctionReturn(0); 2005 } 2006 2007 /* -------------------------------------------------------------- */ 2008 #undef __FUNCT__ 2009 #define __FUNCT__ "MatPermute_SeqAIJ" 2010 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2011 { 2012 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2013 PetscErrorCode ierr; 2014 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2015 const PetscInt *row,*col; 2016 PetscInt *cnew,j,*lens; 2017 IS icolp,irowp; 2018 PetscInt *cwork = PETSC_NULL; 2019 PetscScalar *vwork = PETSC_NULL; 2020 2021 PetscFunctionBegin; 2022 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2023 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2024 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2025 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2026 2027 /* determine lengths of permuted rows */ 2028 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2029 for (i=0; i<m; i++) { 2030 lens[row[i]] = a->i[i+1] - a->i[i]; 2031 } 2032 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2033 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2034 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2035 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2036 ierr = PetscFree(lens);CHKERRQ(ierr); 2037 2038 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2039 for (i=0; i<m; i++) { 2040 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2041 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2042 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2043 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2044 } 2045 ierr = PetscFree(cnew);CHKERRQ(ierr); 2046 (*B)->assembled = PETSC_FALSE; 2047 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2048 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2049 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2050 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2051 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2052 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2053 PetscFunctionReturn(0); 2054 } 2055 2056 #undef __FUNCT__ 2057 #define __FUNCT__ "MatCopy_SeqAIJ" 2058 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2059 { 2060 PetscErrorCode ierr; 2061 2062 PetscFunctionBegin; 2063 /* If the two matrices have the same copy implementation, use fast copy. */ 2064 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2065 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2066 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2067 2068 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2069 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2070 } 2071 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2072 } else { 2073 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2074 } 2075 PetscFunctionReturn(0); 2076 } 2077 2078 #undef __FUNCT__ 2079 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2080 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2081 { 2082 PetscErrorCode ierr; 2083 2084 PetscFunctionBegin; 2085 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2086 PetscFunctionReturn(0); 2087 } 2088 2089 #undef __FUNCT__ 2090 #define __FUNCT__ "MatGetArray_SeqAIJ" 2091 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2092 { 2093 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2094 PetscFunctionBegin; 2095 *array = a->a; 2096 PetscFunctionReturn(0); 2097 } 2098 2099 #undef __FUNCT__ 2100 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2101 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2102 { 2103 PetscFunctionBegin; 2104 PetscFunctionReturn(0); 2105 } 2106 2107 #undef __FUNCT__ 2108 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2109 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2110 { 2111 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2112 PetscErrorCode ierr; 2113 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2114 PetscScalar dx,*y,*xx,*w3_array; 2115 PetscScalar *vscale_array; 2116 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2117 Vec w1,w2,w3; 2118 void *fctx = coloring->fctx; 2119 PetscTruth flg = PETSC_FALSE; 2120 2121 PetscFunctionBegin; 2122 if (!coloring->w1) { 2123 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2124 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2125 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2126 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2127 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2128 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2129 } 2130 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2131 2132 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2133 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2134 if (flg) { 2135 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2136 } else { 2137 PetscTruth assembled; 2138 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2139 if (assembled) { 2140 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2141 } 2142 } 2143 2144 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2145 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2146 2147 /* 2148 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2149 coloring->F for the coarser grids from the finest 2150 */ 2151 if (coloring->F) { 2152 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2153 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2154 if (m1 != m2) { 2155 coloring->F = 0; 2156 } 2157 } 2158 2159 if (coloring->F) { 2160 w1 = coloring->F; 2161 coloring->F = 0; 2162 } else { 2163 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2164 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2165 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2166 } 2167 2168 /* 2169 Compute all the scale factors and share with other processors 2170 */ 2171 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2172 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2173 for (k=0; k<coloring->ncolors; k++) { 2174 /* 2175 Loop over each column associated with color adding the 2176 perturbation to the vector w3. 2177 */ 2178 for (l=0; l<coloring->ncolumns[k]; l++) { 2179 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2180 dx = xx[col]; 2181 if (dx == 0.0) dx = 1.0; 2182 #if !defined(PETSC_USE_COMPLEX) 2183 if (dx < umin && dx >= 0.0) dx = umin; 2184 else if (dx < 0.0 && dx > -umin) dx = -umin; 2185 #else 2186 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2187 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2188 #endif 2189 dx *= epsilon; 2190 vscale_array[col] = 1.0/dx; 2191 } 2192 } 2193 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2194 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2195 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2196 2197 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2198 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2199 2200 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2201 else vscaleforrow = coloring->columnsforrow; 2202 2203 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2204 /* 2205 Loop over each color 2206 */ 2207 for (k=0; k<coloring->ncolors; k++) { 2208 coloring->currentcolor = k; 2209 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2210 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2211 /* 2212 Loop over each column associated with color adding the 2213 perturbation to the vector w3. 2214 */ 2215 for (l=0; l<coloring->ncolumns[k]; l++) { 2216 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2217 dx = xx[col]; 2218 if (dx == 0.0) dx = 1.0; 2219 #if !defined(PETSC_USE_COMPLEX) 2220 if (dx < umin && dx >= 0.0) dx = umin; 2221 else if (dx < 0.0 && dx > -umin) dx = -umin; 2222 #else 2223 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2224 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2225 #endif 2226 dx *= epsilon; 2227 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2228 w3_array[col] += dx; 2229 } 2230 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2231 2232 /* 2233 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2234 */ 2235 2236 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2237 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2238 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2239 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2240 2241 /* 2242 Loop over rows of vector, putting results into Jacobian matrix 2243 */ 2244 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2245 for (l=0; l<coloring->nrows[k]; l++) { 2246 row = coloring->rows[k][l]; 2247 col = coloring->columnsforrow[k][l]; 2248 y[row] *= vscale_array[vscaleforrow[k][l]]; 2249 srow = row + start; 2250 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2251 } 2252 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2253 } 2254 coloring->currentcolor = k; 2255 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2256 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2257 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2258 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2259 PetscFunctionReturn(0); 2260 } 2261 2262 #undef __FUNCT__ 2263 #define __FUNCT__ "MatAXPY_SeqAIJ" 2264 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2265 { 2266 PetscErrorCode ierr; 2267 PetscInt i; 2268 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2269 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2270 2271 PetscFunctionBegin; 2272 if (str == SAME_NONZERO_PATTERN) { 2273 PetscScalar alpha = a; 2274 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2275 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2276 if (y->xtoy && y->XtoY != X) { 2277 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2278 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2279 } 2280 if (!y->xtoy) { /* get xtoy */ 2281 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2282 y->XtoY = X; 2283 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2284 } 2285 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2286 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2287 } else { 2288 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2289 } 2290 PetscFunctionReturn(0); 2291 } 2292 2293 #undef __FUNCT__ 2294 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2295 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2296 { 2297 PetscErrorCode ierr; 2298 2299 PetscFunctionBegin; 2300 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2301 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2302 PetscFunctionReturn(0); 2303 } 2304 2305 #undef __FUNCT__ 2306 #define __FUNCT__ "MatConjugate_SeqAIJ" 2307 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2308 { 2309 #if defined(PETSC_USE_COMPLEX) 2310 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2311 PetscInt i,nz; 2312 PetscScalar *a; 2313 2314 PetscFunctionBegin; 2315 nz = aij->nz; 2316 a = aij->a; 2317 for (i=0; i<nz; i++) { 2318 a[i] = PetscConj(a[i]); 2319 } 2320 #else 2321 PetscFunctionBegin; 2322 #endif 2323 PetscFunctionReturn(0); 2324 } 2325 2326 #undef __FUNCT__ 2327 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2328 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2329 { 2330 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2331 PetscErrorCode ierr; 2332 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2333 PetscReal atmp; 2334 PetscScalar *x; 2335 MatScalar *aa; 2336 2337 PetscFunctionBegin; 2338 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2339 aa = a->a; 2340 ai = a->i; 2341 aj = a->j; 2342 2343 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2344 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2345 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2346 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2347 for (i=0; i<m; i++) { 2348 ncols = ai[1] - ai[0]; ai++; 2349 x[i] = 0.0; 2350 for (j=0; j<ncols; j++){ 2351 atmp = PetscAbsScalar(*aa); 2352 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2353 aa++; aj++; 2354 } 2355 } 2356 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2357 PetscFunctionReturn(0); 2358 } 2359 2360 #undef __FUNCT__ 2361 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2362 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2363 { 2364 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2365 PetscErrorCode ierr; 2366 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2367 PetscScalar *x; 2368 MatScalar *aa; 2369 2370 PetscFunctionBegin; 2371 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2372 aa = a->a; 2373 ai = a->i; 2374 aj = a->j; 2375 2376 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2377 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2378 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2379 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2380 for (i=0; i<m; i++) { 2381 ncols = ai[1] - ai[0]; ai++; 2382 if (ncols == A->cmap->n) { /* row is dense */ 2383 x[i] = *aa; if (idx) idx[i] = 0; 2384 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2385 x[i] = 0.0; 2386 if (idx) { 2387 idx[i] = 0; /* in case ncols is zero */ 2388 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2389 if (aj[j] > j) { 2390 idx[i] = j; 2391 break; 2392 } 2393 } 2394 } 2395 } 2396 for (j=0; j<ncols; j++){ 2397 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2398 aa++; aj++; 2399 } 2400 } 2401 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2402 PetscFunctionReturn(0); 2403 } 2404 2405 #undef __FUNCT__ 2406 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2407 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2408 { 2409 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2410 PetscErrorCode ierr; 2411 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2412 PetscReal atmp; 2413 PetscScalar *x; 2414 MatScalar *aa; 2415 2416 PetscFunctionBegin; 2417 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2418 aa = a->a; 2419 ai = a->i; 2420 aj = a->j; 2421 2422 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2423 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2424 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2425 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2426 for (i=0; i<m; i++) { 2427 ncols = ai[1] - ai[0]; ai++; 2428 if (ncols) { 2429 /* Get first nonzero */ 2430 for(j = 0; j < ncols; j++) { 2431 atmp = PetscAbsScalar(aa[j]); 2432 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2433 } 2434 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2435 } else { 2436 x[i] = 0.0; if (idx) idx[i] = 0; 2437 } 2438 for(j = 0; j < ncols; j++) { 2439 atmp = PetscAbsScalar(*aa); 2440 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2441 aa++; aj++; 2442 } 2443 } 2444 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2445 PetscFunctionReturn(0); 2446 } 2447 2448 #undef __FUNCT__ 2449 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2450 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2451 { 2452 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2453 PetscErrorCode ierr; 2454 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2455 PetscScalar *x; 2456 MatScalar *aa; 2457 2458 PetscFunctionBegin; 2459 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2460 aa = a->a; 2461 ai = a->i; 2462 aj = a->j; 2463 2464 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2465 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2466 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2467 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2468 for (i=0; i<m; i++) { 2469 ncols = ai[1] - ai[0]; ai++; 2470 if (ncols == A->cmap->n) { /* row is dense */ 2471 x[i] = *aa; if (idx) idx[i] = 0; 2472 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2473 x[i] = 0.0; 2474 if (idx) { /* find first implicit 0.0 in the row */ 2475 idx[i] = 0; /* in case ncols is zero */ 2476 for (j=0;j<ncols;j++) { 2477 if (aj[j] > j) { 2478 idx[i] = j; 2479 break; 2480 } 2481 } 2482 } 2483 } 2484 for (j=0; j<ncols; j++){ 2485 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2486 aa++; aj++; 2487 } 2488 } 2489 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2490 PetscFunctionReturn(0); 2491 } 2492 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2493 /* -------------------------------------------------------------------*/ 2494 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2495 MatGetRow_SeqAIJ, 2496 MatRestoreRow_SeqAIJ, 2497 MatMult_SeqAIJ, 2498 /* 4*/ MatMultAdd_SeqAIJ, 2499 MatMultTranspose_SeqAIJ, 2500 MatMultTransposeAdd_SeqAIJ, 2501 0, 2502 0, 2503 0, 2504 /*10*/ 0, 2505 MatLUFactor_SeqAIJ, 2506 0, 2507 MatSOR_SeqAIJ, 2508 MatTranspose_SeqAIJ, 2509 /*15*/ MatGetInfo_SeqAIJ, 2510 MatEqual_SeqAIJ, 2511 MatGetDiagonal_SeqAIJ, 2512 MatDiagonalScale_SeqAIJ, 2513 MatNorm_SeqAIJ, 2514 /*20*/ 0, 2515 MatAssemblyEnd_SeqAIJ, 2516 MatSetOption_SeqAIJ, 2517 MatZeroEntries_SeqAIJ, 2518 /*24*/ MatZeroRows_SeqAIJ, 2519 0, 2520 0, 2521 0, 2522 0, 2523 /*29*/ MatSetUpPreallocation_SeqAIJ, 2524 0, 2525 0, 2526 MatGetArray_SeqAIJ, 2527 MatRestoreArray_SeqAIJ, 2528 /*34*/ MatDuplicate_SeqAIJ, 2529 0, 2530 0, 2531 MatILUFactor_SeqAIJ, 2532 0, 2533 /*39*/ MatAXPY_SeqAIJ, 2534 MatGetSubMatrices_SeqAIJ, 2535 MatIncreaseOverlap_SeqAIJ, 2536 MatGetValues_SeqAIJ, 2537 MatCopy_SeqAIJ, 2538 /*44*/ MatGetRowMax_SeqAIJ, 2539 MatScale_SeqAIJ, 2540 0, 2541 MatDiagonalSet_SeqAIJ, 2542 0, 2543 /*49*/ MatSetBlockSize_SeqAIJ, 2544 MatGetRowIJ_SeqAIJ, 2545 MatRestoreRowIJ_SeqAIJ, 2546 MatGetColumnIJ_SeqAIJ, 2547 MatRestoreColumnIJ_SeqAIJ, 2548 /*54*/ MatFDColoringCreate_SeqAIJ, 2549 0, 2550 0, 2551 MatPermute_SeqAIJ, 2552 0, 2553 /*59*/ 0, 2554 MatDestroy_SeqAIJ, 2555 MatView_SeqAIJ, 2556 0, 2557 0, 2558 /*64*/ 0, 2559 0, 2560 0, 2561 0, 2562 0, 2563 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2564 MatGetRowMinAbs_SeqAIJ, 2565 0, 2566 MatSetColoring_SeqAIJ, 2567 #if defined(PETSC_HAVE_ADIC) 2568 MatSetValuesAdic_SeqAIJ, 2569 #else 2570 0, 2571 #endif 2572 /*74*/ MatSetValuesAdifor_SeqAIJ, 2573 MatFDColoringApply_AIJ, 2574 0, 2575 0, 2576 0, 2577 /*79*/ 0, 2578 0, 2579 0, 2580 0, 2581 MatLoad_SeqAIJ, 2582 /*84*/ MatIsSymmetric_SeqAIJ, 2583 MatIsHermitian_SeqAIJ, 2584 0, 2585 0, 2586 0, 2587 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2588 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2589 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2590 MatPtAP_Basic, 2591 MatPtAPSymbolic_SeqAIJ, 2592 /*94*/ MatPtAPNumeric_SeqAIJ, 2593 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2594 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2595 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2596 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2597 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2598 0, 2599 0, 2600 MatConjugate_SeqAIJ, 2601 0, 2602 /*104*/MatSetValuesRow_SeqAIJ, 2603 MatRealPart_SeqAIJ, 2604 MatImaginaryPart_SeqAIJ, 2605 0, 2606 0, 2607 /*109*/0, 2608 0, 2609 MatGetRowMin_SeqAIJ, 2610 0, 2611 MatMissingDiagonal_SeqAIJ, 2612 /*114*/0, 2613 0, 2614 0, 2615 0, 2616 0 2617 }; 2618 2619 EXTERN_C_BEGIN 2620 #undef __FUNCT__ 2621 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2622 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2623 { 2624 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2625 PetscInt i,nz,n; 2626 2627 PetscFunctionBegin; 2628 2629 nz = aij->maxnz; 2630 n = mat->rmap->n; 2631 for (i=0; i<nz; i++) { 2632 aij->j[i] = indices[i]; 2633 } 2634 aij->nz = nz; 2635 for (i=0; i<n; i++) { 2636 aij->ilen[i] = aij->imax[i]; 2637 } 2638 2639 PetscFunctionReturn(0); 2640 } 2641 EXTERN_C_END 2642 2643 #undef __FUNCT__ 2644 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2645 /*@ 2646 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2647 in the matrix. 2648 2649 Input Parameters: 2650 + mat - the SeqAIJ matrix 2651 - indices - the column indices 2652 2653 Level: advanced 2654 2655 Notes: 2656 This can be called if you have precomputed the nonzero structure of the 2657 matrix and want to provide it to the matrix object to improve the performance 2658 of the MatSetValues() operation. 2659 2660 You MUST have set the correct numbers of nonzeros per row in the call to 2661 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2662 2663 MUST be called before any calls to MatSetValues(); 2664 2665 The indices should start with zero, not one. 2666 2667 @*/ 2668 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2669 { 2670 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2671 2672 PetscFunctionBegin; 2673 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2674 PetscValidPointer(indices,2); 2675 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2676 if (f) { 2677 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2678 } else { 2679 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2680 } 2681 PetscFunctionReturn(0); 2682 } 2683 2684 /* ----------------------------------------------------------------------------------------*/ 2685 2686 EXTERN_C_BEGIN 2687 #undef __FUNCT__ 2688 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2689 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2690 { 2691 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2692 PetscErrorCode ierr; 2693 size_t nz = aij->i[mat->rmap->n]; 2694 2695 PetscFunctionBegin; 2696 if (aij->nonew != 1) { 2697 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2698 } 2699 2700 /* allocate space for values if not already there */ 2701 if (!aij->saved_values) { 2702 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2703 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2704 } 2705 2706 /* copy values over */ 2707 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2708 PetscFunctionReturn(0); 2709 } 2710 EXTERN_C_END 2711 2712 #undef __FUNCT__ 2713 #define __FUNCT__ "MatStoreValues" 2714 /*@ 2715 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2716 example, reuse of the linear part of a Jacobian, while recomputing the 2717 nonlinear portion. 2718 2719 Collect on Mat 2720 2721 Input Parameters: 2722 . mat - the matrix (currently only AIJ matrices support this option) 2723 2724 Level: advanced 2725 2726 Common Usage, with SNESSolve(): 2727 $ Create Jacobian matrix 2728 $ Set linear terms into matrix 2729 $ Apply boundary conditions to matrix, at this time matrix must have 2730 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2731 $ boundary conditions again will not change the nonzero structure 2732 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2733 $ ierr = MatStoreValues(mat); 2734 $ Call SNESSetJacobian() with matrix 2735 $ In your Jacobian routine 2736 $ ierr = MatRetrieveValues(mat); 2737 $ Set nonlinear terms in matrix 2738 2739 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2740 $ // build linear portion of Jacobian 2741 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2742 $ ierr = MatStoreValues(mat); 2743 $ loop over nonlinear iterations 2744 $ ierr = MatRetrieveValues(mat); 2745 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2746 $ // call MatAssemblyBegin/End() on matrix 2747 $ Solve linear system with Jacobian 2748 $ endloop 2749 2750 Notes: 2751 Matrix must already be assemblied before calling this routine 2752 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2753 calling this routine. 2754 2755 When this is called multiple times it overwrites the previous set of stored values 2756 and does not allocated additional space. 2757 2758 .seealso: MatRetrieveValues() 2759 2760 @*/ 2761 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2762 { 2763 PetscErrorCode ierr,(*f)(Mat); 2764 2765 PetscFunctionBegin; 2766 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2767 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2768 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2769 2770 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2771 if (f) { 2772 ierr = (*f)(mat);CHKERRQ(ierr); 2773 } else { 2774 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2775 } 2776 PetscFunctionReturn(0); 2777 } 2778 2779 EXTERN_C_BEGIN 2780 #undef __FUNCT__ 2781 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2782 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2783 { 2784 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2785 PetscErrorCode ierr; 2786 PetscInt nz = aij->i[mat->rmap->n]; 2787 2788 PetscFunctionBegin; 2789 if (aij->nonew != 1) { 2790 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2791 } 2792 if (!aij->saved_values) { 2793 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2794 } 2795 /* copy values over */ 2796 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2797 PetscFunctionReturn(0); 2798 } 2799 EXTERN_C_END 2800 2801 #undef __FUNCT__ 2802 #define __FUNCT__ "MatRetrieveValues" 2803 /*@ 2804 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2805 example, reuse of the linear part of a Jacobian, while recomputing the 2806 nonlinear portion. 2807 2808 Collect on Mat 2809 2810 Input Parameters: 2811 . mat - the matrix (currently on AIJ matrices support this option) 2812 2813 Level: advanced 2814 2815 .seealso: MatStoreValues() 2816 2817 @*/ 2818 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2819 { 2820 PetscErrorCode ierr,(*f)(Mat); 2821 2822 PetscFunctionBegin; 2823 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2824 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2825 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2826 2827 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2828 if (f) { 2829 ierr = (*f)(mat);CHKERRQ(ierr); 2830 } else { 2831 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2832 } 2833 PetscFunctionReturn(0); 2834 } 2835 2836 2837 /* --------------------------------------------------------------------------------*/ 2838 #undef __FUNCT__ 2839 #define __FUNCT__ "MatCreateSeqAIJ" 2840 /*@C 2841 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2842 (the default parallel PETSc format). For good matrix assembly performance 2843 the user should preallocate the matrix storage by setting the parameter nz 2844 (or the array nnz). By setting these parameters accurately, performance 2845 during matrix assembly can be increased by more than a factor of 50. 2846 2847 Collective on MPI_Comm 2848 2849 Input Parameters: 2850 + comm - MPI communicator, set to PETSC_COMM_SELF 2851 . m - number of rows 2852 . n - number of columns 2853 . nz - number of nonzeros per row (same for all rows) 2854 - nnz - array containing the number of nonzeros in the various rows 2855 (possibly different for each row) or PETSC_NULL 2856 2857 Output Parameter: 2858 . A - the matrix 2859 2860 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2861 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2862 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2863 2864 Notes: 2865 If nnz is given then nz is ignored 2866 2867 The AIJ format (also called the Yale sparse matrix format or 2868 compressed row storage), is fully compatible with standard Fortran 77 2869 storage. That is, the stored row and column indices can begin at 2870 either one (as in Fortran) or zero. See the users' manual for details. 2871 2872 Specify the preallocated storage with either nz or nnz (not both). 2873 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2874 allocation. For large problems you MUST preallocate memory or you 2875 will get TERRIBLE performance, see the users' manual chapter on matrices. 2876 2877 By default, this format uses inodes (identical nodes) when possible, to 2878 improve numerical efficiency of matrix-vector products and solves. We 2879 search for consecutive rows with the same nonzero structure, thereby 2880 reusing matrix information to achieve increased efficiency. 2881 2882 Options Database Keys: 2883 + -mat_no_inode - Do not use inodes 2884 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2885 - -mat_aij_oneindex - Internally use indexing starting at 1 2886 rather than 0. Note that when calling MatSetValues(), 2887 the user still MUST index entries starting at 0! 2888 2889 Level: intermediate 2890 2891 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2892 2893 @*/ 2894 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2895 { 2896 PetscErrorCode ierr; 2897 2898 PetscFunctionBegin; 2899 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2900 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2901 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2902 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2903 PetscFunctionReturn(0); 2904 } 2905 2906 #undef __FUNCT__ 2907 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2908 /*@C 2909 MatSeqAIJSetPreallocation - For good matrix assembly performance 2910 the user should preallocate the matrix storage by setting the parameter nz 2911 (or the array nnz). By setting these parameters accurately, performance 2912 during matrix assembly can be increased by more than a factor of 50. 2913 2914 Collective on MPI_Comm 2915 2916 Input Parameters: 2917 + B - The matrix-free 2918 . nz - number of nonzeros per row (same for all rows) 2919 - nnz - array containing the number of nonzeros in the various rows 2920 (possibly different for each row) or PETSC_NULL 2921 2922 Notes: 2923 If nnz is given then nz is ignored 2924 2925 The AIJ format (also called the Yale sparse matrix format or 2926 compressed row storage), is fully compatible with standard Fortran 77 2927 storage. That is, the stored row and column indices can begin at 2928 either one (as in Fortran) or zero. See the users' manual for details. 2929 2930 Specify the preallocated storage with either nz or nnz (not both). 2931 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2932 allocation. For large problems you MUST preallocate memory or you 2933 will get TERRIBLE performance, see the users' manual chapter on matrices. 2934 2935 You can call MatGetInfo() to get information on how effective the preallocation was; 2936 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2937 You can also run with the option -info and look for messages with the string 2938 malloc in them to see if additional memory allocation was needed. 2939 2940 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2941 entries or columns indices 2942 2943 By default, this format uses inodes (identical nodes) when possible, to 2944 improve numerical efficiency of matrix-vector products and solves. We 2945 search for consecutive rows with the same nonzero structure, thereby 2946 reusing matrix information to achieve increased efficiency. 2947 2948 Options Database Keys: 2949 + -mat_no_inode - Do not use inodes 2950 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2951 - -mat_aij_oneindex - Internally use indexing starting at 1 2952 rather than 0. Note that when calling MatSetValues(), 2953 the user still MUST index entries starting at 0! 2954 2955 Level: intermediate 2956 2957 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2958 2959 @*/ 2960 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2961 { 2962 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2963 2964 PetscFunctionBegin; 2965 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2966 if (f) { 2967 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2968 } 2969 PetscFunctionReturn(0); 2970 } 2971 2972 EXTERN_C_BEGIN 2973 #undef __FUNCT__ 2974 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2975 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2976 { 2977 Mat_SeqAIJ *b; 2978 PetscTruth skipallocation = PETSC_FALSE; 2979 PetscErrorCode ierr; 2980 PetscInt i; 2981 2982 PetscFunctionBegin; 2983 2984 if (nz == MAT_SKIP_ALLOCATION) { 2985 skipallocation = PETSC_TRUE; 2986 nz = 0; 2987 } 2988 2989 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 2990 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 2991 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2992 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2993 2994 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2995 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2996 if (nnz) { 2997 for (i=0; i<B->rmap->n; i++) { 2998 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2999 if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap->n); 3000 } 3001 } 3002 3003 B->preallocated = PETSC_TRUE; 3004 b = (Mat_SeqAIJ*)B->data; 3005 3006 if (!skipallocation) { 3007 if (!b->imax) { 3008 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3009 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3010 } 3011 if (!nnz) { 3012 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3013 else if (nz <= 0) nz = 1; 3014 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3015 nz = nz*B->rmap->n; 3016 } else { 3017 nz = 0; 3018 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3019 } 3020 /* b->ilen will count nonzeros in each row so far. */ 3021 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3022 3023 /* allocate the matrix space */ 3024 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3025 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3026 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3027 b->i[0] = 0; 3028 for (i=1; i<B->rmap->n+1; i++) { 3029 b->i[i] = b->i[i-1] + b->imax[i-1]; 3030 } 3031 b->singlemalloc = PETSC_TRUE; 3032 b->free_a = PETSC_TRUE; 3033 b->free_ij = PETSC_TRUE; 3034 } else { 3035 b->free_a = PETSC_FALSE; 3036 b->free_ij = PETSC_FALSE; 3037 } 3038 3039 b->nz = 0; 3040 b->maxnz = nz; 3041 B->info.nz_unneeded = (double)b->maxnz; 3042 PetscFunctionReturn(0); 3043 } 3044 EXTERN_C_END 3045 3046 #undef __FUNCT__ 3047 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3048 /*@ 3049 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3050 3051 Input Parameters: 3052 + B - the matrix 3053 . i - the indices into j for the start of each row (starts with zero) 3054 . j - the column indices for each row (starts with zero) these must be sorted for each row 3055 - v - optional values in the matrix 3056 3057 Level: developer 3058 3059 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3060 3061 .keywords: matrix, aij, compressed row, sparse, sequential 3062 3063 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3064 @*/ 3065 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3066 { 3067 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3068 PetscErrorCode ierr; 3069 3070 PetscFunctionBegin; 3071 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3072 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3073 if (f) { 3074 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3075 } 3076 PetscFunctionReturn(0); 3077 } 3078 3079 EXTERN_C_BEGIN 3080 #undef __FUNCT__ 3081 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3082 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3083 { 3084 PetscInt i; 3085 PetscInt m,n; 3086 PetscInt nz; 3087 PetscInt *nnz, nz_max = 0; 3088 PetscScalar *values; 3089 PetscErrorCode ierr; 3090 3091 PetscFunctionBegin; 3092 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3093 3094 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3095 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3096 for(i = 0; i < m; i++) { 3097 nz = Ii[i+1]- Ii[i]; 3098 nz_max = PetscMax(nz_max, nz); 3099 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3100 nnz[i] = nz; 3101 } 3102 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3103 ierr = PetscFree(nnz);CHKERRQ(ierr); 3104 3105 if (v) { 3106 values = (PetscScalar*) v; 3107 } else { 3108 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3109 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3110 } 3111 3112 for(i = 0; i < m; i++) { 3113 nz = Ii[i+1] - Ii[i]; 3114 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3115 } 3116 3117 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3118 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3119 3120 if (!v) { 3121 ierr = PetscFree(values);CHKERRQ(ierr); 3122 } 3123 PetscFunctionReturn(0); 3124 } 3125 EXTERN_C_END 3126 3127 #include "../src/mat/impls/dense/seq/dense.h" 3128 #include "private/petscaxpy.h" 3129 3130 #undef __FUNCT__ 3131 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3132 /* 3133 Computes (B'*A')' since computing B*A directly is untenable 3134 3135 n p p 3136 ( ) ( ) ( ) 3137 m ( A ) * n ( B ) = m ( C ) 3138 ( ) ( ) ( ) 3139 3140 */ 3141 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3142 { 3143 PetscErrorCode ierr; 3144 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3145 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3146 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3147 PetscInt i,n,m,q,p; 3148 const PetscInt *ii,*idx; 3149 const PetscScalar *b,*a,*a_q; 3150 PetscScalar *c,*c_q; 3151 3152 PetscFunctionBegin; 3153 m = A->rmap->n; 3154 n = A->cmap->n; 3155 p = B->cmap->n; 3156 a = sub_a->v; 3157 b = sub_b->a; 3158 c = sub_c->v; 3159 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3160 3161 ii = sub_b->i; 3162 idx = sub_b->j; 3163 for (i=0; i<n; i++) { 3164 q = ii[i+1] - ii[i]; 3165 while (q-->0) { 3166 c_q = c + m*(*idx); 3167 a_q = a + m*i; 3168 PetscAXPY(c_q,*b,a_q,m); 3169 idx++; 3170 b++; 3171 } 3172 } 3173 PetscFunctionReturn(0); 3174 } 3175 3176 #undef __FUNCT__ 3177 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3178 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3179 { 3180 PetscErrorCode ierr; 3181 PetscInt m=A->rmap->n,n=B->cmap->n; 3182 Mat Cmat; 3183 3184 PetscFunctionBegin; 3185 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 3186 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3187 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3188 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3189 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3190 Cmat->assembled = PETSC_TRUE; 3191 *C = Cmat; 3192 PetscFunctionReturn(0); 3193 } 3194 3195 /* ----------------------------------------------------------------*/ 3196 #undef __FUNCT__ 3197 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3198 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3199 { 3200 PetscErrorCode ierr; 3201 3202 PetscFunctionBegin; 3203 if (scall == MAT_INITIAL_MATRIX){ 3204 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3205 } 3206 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3207 PetscFunctionReturn(0); 3208 } 3209 3210 3211 /*MC 3212 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3213 based on compressed sparse row format. 3214 3215 Options Database Keys: 3216 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3217 3218 Level: beginner 3219 3220 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3221 M*/ 3222 3223 EXTERN_C_BEGIN 3224 #if defined(PETSC_HAVE_PASTIX) 3225 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3226 #endif 3227 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3228 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3229 #endif 3230 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3231 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3232 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3233 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3234 #if defined(PETSC_HAVE_MUMPS) 3235 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3236 #endif 3237 #if defined(PETSC_HAVE_SUPERLU) 3238 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3239 #endif 3240 #if defined(PETSC_HAVE_SUPERLU_DIST) 3241 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3242 #endif 3243 #if defined(PETSC_HAVE_SPOOLES) 3244 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3245 #endif 3246 #if defined(PETSC_HAVE_UMFPACK) 3247 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3248 #endif 3249 #if defined(PETSC_HAVE_CHOLMOD) 3250 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3251 #endif 3252 #if defined(PETSC_HAVE_LUSOL) 3253 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3254 #endif 3255 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3256 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3257 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3258 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3259 #endif 3260 EXTERN_C_END 3261 3262 3263 EXTERN_C_BEGIN 3264 #undef __FUNCT__ 3265 #define __FUNCT__ "MatCreate_SeqAIJ" 3266 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3267 { 3268 Mat_SeqAIJ *b; 3269 PetscErrorCode ierr; 3270 PetscMPIInt size; 3271 3272 PetscFunctionBegin; 3273 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3274 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3275 3276 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3277 B->data = (void*)b; 3278 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3279 B->mapping = 0; 3280 b->row = 0; 3281 b->col = 0; 3282 b->icol = 0; 3283 b->reallocs = 0; 3284 b->ignorezeroentries = PETSC_FALSE; 3285 b->roworiented = PETSC_TRUE; 3286 b->nonew = 0; 3287 b->diag = 0; 3288 b->solve_work = 0; 3289 B->spptr = 0; 3290 b->saved_values = 0; 3291 b->idiag = 0; 3292 b->mdiag = 0; 3293 b->ssor_work = 0; 3294 b->omega = 1.0; 3295 b->fshift = 0.0; 3296 b->idiagvalid = PETSC_FALSE; 3297 b->keepnonzeropattern = PETSC_FALSE; 3298 b->xtoy = 0; 3299 b->XtoY = 0; 3300 b->compressedrow.use = PETSC_FALSE; 3301 b->compressedrow.nrows = B->rmap->n; 3302 b->compressedrow.i = PETSC_NULL; 3303 b->compressedrow.rindex = PETSC_NULL; 3304 b->compressedrow.checked = PETSC_FALSE; 3305 B->same_nonzero = PETSC_FALSE; 3306 3307 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3308 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3309 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3310 "MatGetFactor_seqaij_matlab", 3311 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3312 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3313 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3314 #endif 3315 #if defined(PETSC_HAVE_PASTIX) 3316 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3317 "MatGetFactor_seqaij_pastix", 3318 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3319 #endif 3320 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3321 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3322 "MatGetFactor_seqaij_essl", 3323 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3324 #endif 3325 #if defined(PETSC_HAVE_SUPERLU) 3326 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3327 "MatGetFactor_seqaij_superlu", 3328 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3329 #endif 3330 #if defined(PETSC_HAVE_SUPERLU_DIST) 3331 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3332 "MatGetFactor_seqaij_superlu_dist", 3333 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3334 #endif 3335 #if defined(PETSC_HAVE_SPOOLES) 3336 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3337 "MatGetFactor_seqaij_spooles", 3338 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3339 #endif 3340 #if defined(PETSC_HAVE_MUMPS) 3341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3342 "MatGetFactor_aij_mumps", 3343 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3344 #endif 3345 #if defined(PETSC_HAVE_UMFPACK) 3346 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3347 "MatGetFactor_seqaij_umfpack", 3348 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3349 #endif 3350 #if defined(PETSC_HAVE_CHOLMOD) 3351 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3352 "MatGetFactor_seqaij_cholmod", 3353 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3354 #endif 3355 #if defined(PETSC_HAVE_LUSOL) 3356 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3357 "MatGetFactor_seqaij_lusol", 3358 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3359 #endif 3360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3361 "MatGetFactor_seqaij_petsc", 3362 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3363 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3364 "MatGetFactorAvailable_seqaij_petsc", 3365 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3366 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3367 "MatGetFactor_seqaij_bas", 3368 MatGetFactor_seqaij_bas); 3369 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3370 "MatSeqAIJSetColumnIndices_SeqAIJ", 3371 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3372 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3373 "MatStoreValues_SeqAIJ", 3374 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3375 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3376 "MatRetrieveValues_SeqAIJ", 3377 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3378 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3379 "MatConvert_SeqAIJ_SeqSBAIJ", 3380 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3381 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3382 "MatConvert_SeqAIJ_SeqBAIJ", 3383 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3384 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3385 "MatConvert_SeqAIJ_SeqCSRPERM", 3386 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3387 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3388 "MatConvert_SeqAIJ_SeqCRL", 3389 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3390 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3391 "MatIsTranspose_SeqAIJ", 3392 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3393 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3394 "MatIsHermitianTranspose_SeqAIJ", 3395 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3396 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3397 "MatSeqAIJSetPreallocation_SeqAIJ", 3398 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3399 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3400 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3401 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3402 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3403 "MatReorderForNonzeroDiagonal_SeqAIJ", 3404 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3405 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3406 "MatMatMult_SeqDense_SeqAIJ", 3407 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3408 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3409 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3410 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3411 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3412 "MatMatMultNumeric_SeqDense_SeqAIJ", 3413 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3414 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3415 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3416 PetscFunctionReturn(0); 3417 } 3418 EXTERN_C_END 3419 3420 #undef __FUNCT__ 3421 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3422 /* 3423 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3424 */ 3425 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3426 { 3427 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3428 PetscErrorCode ierr; 3429 PetscInt i,m = A->rmap->n; 3430 3431 PetscFunctionBegin; 3432 c = (Mat_SeqAIJ*)C->data; 3433 3434 C->factortype = A->factortype; 3435 c->row = 0; 3436 c->col = 0; 3437 c->icol = 0; 3438 c->reallocs = 0; 3439 3440 C->assembled = PETSC_TRUE; 3441 3442 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3443 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3444 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3445 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3446 3447 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3448 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3449 for (i=0; i<m; i++) { 3450 c->imax[i] = a->imax[i]; 3451 c->ilen[i] = a->ilen[i]; 3452 } 3453 3454 /* allocate the matrix space */ 3455 if (mallocmatspace){ 3456 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3457 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3458 c->singlemalloc = PETSC_TRUE; 3459 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3460 if (m > 0) { 3461 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3462 if (cpvalues == MAT_COPY_VALUES) { 3463 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3464 } else { 3465 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3466 } 3467 } 3468 } 3469 3470 c->ignorezeroentries = a->ignorezeroentries; 3471 c->roworiented = a->roworiented; 3472 c->nonew = a->nonew; 3473 if (a->diag) { 3474 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3475 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3476 for (i=0; i<m; i++) { 3477 c->diag[i] = a->diag[i]; 3478 } 3479 } else c->diag = 0; 3480 c->solve_work = 0; 3481 c->saved_values = 0; 3482 c->idiag = 0; 3483 c->ssor_work = 0; 3484 c->keepnonzeropattern = a->keepnonzeropattern; 3485 c->free_a = PETSC_TRUE; 3486 c->free_ij = PETSC_TRUE; 3487 c->xtoy = 0; 3488 c->XtoY = 0; 3489 3490 c->nz = a->nz; 3491 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3492 C->preallocated = PETSC_TRUE; 3493 3494 c->compressedrow.use = a->compressedrow.use; 3495 c->compressedrow.nrows = a->compressedrow.nrows; 3496 c->compressedrow.checked = a->compressedrow.checked; 3497 if (a->compressedrow.checked && a->compressedrow.use){ 3498 i = a->compressedrow.nrows; 3499 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3500 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3501 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3502 } else { 3503 c->compressedrow.use = PETSC_FALSE; 3504 c->compressedrow.i = PETSC_NULL; 3505 c->compressedrow.rindex = PETSC_NULL; 3506 } 3507 C->same_nonzero = A->same_nonzero; 3508 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3509 3510 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3511 PetscFunctionReturn(0); 3512 } 3513 3514 #undef __FUNCT__ 3515 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3516 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3517 { 3518 PetscErrorCode ierr; 3519 3520 PetscFunctionBegin; 3521 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3522 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3523 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3524 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3525 PetscFunctionReturn(0); 3526 } 3527 3528 #undef __FUNCT__ 3529 #define __FUNCT__ "MatLoad_SeqAIJ" 3530 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3531 { 3532 Mat_SeqAIJ *a; 3533 Mat B; 3534 PetscErrorCode ierr; 3535 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3536 int fd; 3537 PetscMPIInt size; 3538 MPI_Comm comm; 3539 3540 PetscFunctionBegin; 3541 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3542 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3543 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3544 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3545 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3546 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3547 M = header[1]; N = header[2]; nz = header[3]; 3548 3549 if (nz < 0) { 3550 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3551 } 3552 3553 /* read in row lengths */ 3554 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3555 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3556 3557 /* check if sum of rowlengths is same as nz */ 3558 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3559 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3560 3561 /* create our matrix */ 3562 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3563 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3564 ierr = MatSetType(B,type);CHKERRQ(ierr); 3565 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3566 a = (Mat_SeqAIJ*)B->data; 3567 3568 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3569 3570 /* read in nonzero values */ 3571 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3572 3573 /* set matrix "i" values */ 3574 a->i[0] = 0; 3575 for (i=1; i<= M; i++) { 3576 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3577 a->ilen[i-1] = rowlengths[i-1]; 3578 } 3579 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3580 3581 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3582 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3583 *A = B; 3584 PetscFunctionReturn(0); 3585 } 3586 3587 #undef __FUNCT__ 3588 #define __FUNCT__ "MatEqual_SeqAIJ" 3589 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3590 { 3591 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3592 PetscErrorCode ierr; 3593 #if defined(PETSC_USE_COMPLEX) 3594 PetscInt k; 3595 #endif 3596 3597 PetscFunctionBegin; 3598 /* If the matrix dimensions are not equal,or no of nonzeros */ 3599 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3600 *flg = PETSC_FALSE; 3601 PetscFunctionReturn(0); 3602 } 3603 3604 /* if the a->i are the same */ 3605 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3606 if (!*flg) PetscFunctionReturn(0); 3607 3608 /* if a->j are the same */ 3609 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3610 if (!*flg) PetscFunctionReturn(0); 3611 3612 /* if a->a are the same */ 3613 #if defined(PETSC_USE_COMPLEX) 3614 for (k=0; k<a->nz; k++){ 3615 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3616 *flg = PETSC_FALSE; 3617 PetscFunctionReturn(0); 3618 } 3619 } 3620 #else 3621 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3622 #endif 3623 PetscFunctionReturn(0); 3624 } 3625 3626 #undef __FUNCT__ 3627 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3628 /*@ 3629 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3630 provided by the user. 3631 3632 Collective on MPI_Comm 3633 3634 Input Parameters: 3635 + comm - must be an MPI communicator of size 1 3636 . m - number of rows 3637 . n - number of columns 3638 . i - row indices 3639 . j - column indices 3640 - a - matrix values 3641 3642 Output Parameter: 3643 . mat - the matrix 3644 3645 Level: intermediate 3646 3647 Notes: 3648 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3649 once the matrix is destroyed 3650 3651 You cannot set new nonzero locations into this matrix, that will generate an error. 3652 3653 The i and j indices are 0 based 3654 3655 The format which is used for the sparse matrix input, is equivalent to a 3656 row-major ordering.. i.e for the following matrix, the input data expected is 3657 as shown: 3658 3659 1 0 0 3660 2 0 3 3661 4 5 6 3662 3663 i = {0,1,3,6} [size = nrow+1 = 3+1] 3664 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3665 v = {1,2,3,4,5,6} [size = nz = 6] 3666 3667 3668 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3669 3670 @*/ 3671 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3672 { 3673 PetscErrorCode ierr; 3674 PetscInt ii; 3675 Mat_SeqAIJ *aij; 3676 #if defined(PETSC_USE_DEBUG) 3677 PetscInt jj; 3678 #endif 3679 3680 PetscFunctionBegin; 3681 if (i[0]) { 3682 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3683 } 3684 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3685 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3686 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3687 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3688 aij = (Mat_SeqAIJ*)(*mat)->data; 3689 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3690 3691 aij->i = i; 3692 aij->j = j; 3693 aij->a = a; 3694 aij->singlemalloc = PETSC_FALSE; 3695 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3696 aij->free_a = PETSC_FALSE; 3697 aij->free_ij = PETSC_FALSE; 3698 3699 for (ii=0; ii<m; ii++) { 3700 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3701 #if defined(PETSC_USE_DEBUG) 3702 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 3703 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3704 if (j[jj] < j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is not sorted",jj-i[ii],j[jj],ii); 3705 if (j[jj] == j[jj]-1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii); 3706 } 3707 #endif 3708 } 3709 #if defined(PETSC_USE_DEBUG) 3710 for (ii=0; ii<aij->i[m]; ii++) { 3711 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3712 if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3713 } 3714 #endif 3715 3716 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3717 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3718 PetscFunctionReturn(0); 3719 } 3720 3721 #undef __FUNCT__ 3722 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3723 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3724 { 3725 PetscErrorCode ierr; 3726 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3727 3728 PetscFunctionBegin; 3729 if (coloring->ctype == IS_COLORING_GLOBAL) { 3730 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3731 a->coloring = coloring; 3732 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3733 PetscInt i,*larray; 3734 ISColoring ocoloring; 3735 ISColoringValue *colors; 3736 3737 /* set coloring for diagonal portion */ 3738 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3739 for (i=0; i<A->cmap->n; i++) { 3740 larray[i] = i; 3741 } 3742 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3743 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3744 for (i=0; i<A->cmap->n; i++) { 3745 colors[i] = coloring->colors[larray[i]]; 3746 } 3747 ierr = PetscFree(larray);CHKERRQ(ierr); 3748 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3749 a->coloring = ocoloring; 3750 } 3751 PetscFunctionReturn(0); 3752 } 3753 3754 #if defined(PETSC_HAVE_ADIC) 3755 EXTERN_C_BEGIN 3756 #include "adic/ad_utils.h" 3757 EXTERN_C_END 3758 3759 #undef __FUNCT__ 3760 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3761 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3762 { 3763 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3764 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3765 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3766 ISColoringValue *color; 3767 3768 PetscFunctionBegin; 3769 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3770 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3771 color = a->coloring->colors; 3772 /* loop over rows */ 3773 for (i=0; i<m; i++) { 3774 nz = ii[i+1] - ii[i]; 3775 /* loop over columns putting computed value into matrix */ 3776 for (j=0; j<nz; j++) { 3777 *v++ = values[color[*jj++]]; 3778 } 3779 values += nlen; /* jump to next row of derivatives */ 3780 } 3781 PetscFunctionReturn(0); 3782 } 3783 #endif 3784 3785 #undef __FUNCT__ 3786 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3787 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3788 { 3789 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3790 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3791 MatScalar *v = a->a; 3792 PetscScalar *values = (PetscScalar *)advalues; 3793 ISColoringValue *color; 3794 3795 PetscFunctionBegin; 3796 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3797 color = a->coloring->colors; 3798 /* loop over rows */ 3799 for (i=0; i<m; i++) { 3800 nz = ii[i+1] - ii[i]; 3801 /* loop over columns putting computed value into matrix */ 3802 for (j=0; j<nz; j++) { 3803 *v++ = values[color[*jj++]]; 3804 } 3805 values += nl; /* jump to next row of derivatives */ 3806 } 3807 PetscFunctionReturn(0); 3808 } 3809 3810 /* 3811 Special version for direct calls from Fortran 3812 */ 3813 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3814 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3815 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3816 #define matsetvaluesseqaij_ matsetvaluesseqaij 3817 #endif 3818 3819 /* Change these macros so can be used in void function */ 3820 #undef CHKERRQ 3821 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3822 #undef SETERRQ2 3823 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3824 3825 EXTERN_C_BEGIN 3826 #undef __FUNCT__ 3827 #define __FUNCT__ "matsetvaluesseqaij_" 3828 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3829 { 3830 Mat A = *AA; 3831 PetscInt m = *mm, n = *nn; 3832 InsertMode is = *isis; 3833 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3834 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3835 PetscInt *imax,*ai,*ailen; 3836 PetscErrorCode ierr; 3837 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3838 MatScalar *ap,value,*aa; 3839 PetscTruth ignorezeroentries = a->ignorezeroentries; 3840 PetscTruth roworiented = a->roworiented; 3841 3842 PetscFunctionBegin; 3843 ierr = MatPreallocated(A);CHKERRQ(ierr); 3844 imax = a->imax; 3845 ai = a->i; 3846 ailen = a->ilen; 3847 aj = a->j; 3848 aa = a->a; 3849 3850 for (k=0; k<m; k++) { /* loop over added rows */ 3851 row = im[k]; 3852 if (row < 0) continue; 3853 #if defined(PETSC_USE_DEBUG) 3854 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3855 #endif 3856 rp = aj + ai[row]; ap = aa + ai[row]; 3857 rmax = imax[row]; nrow = ailen[row]; 3858 low = 0; 3859 high = nrow; 3860 for (l=0; l<n; l++) { /* loop over added columns */ 3861 if (in[l] < 0) continue; 3862 #if defined(PETSC_USE_DEBUG) 3863 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3864 #endif 3865 col = in[l]; 3866 if (roworiented) { 3867 value = v[l + k*n]; 3868 } else { 3869 value = v[k + l*m]; 3870 } 3871 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3872 3873 if (col <= lastcol) low = 0; else high = nrow; 3874 lastcol = col; 3875 while (high-low > 5) { 3876 t = (low+high)/2; 3877 if (rp[t] > col) high = t; 3878 else low = t; 3879 } 3880 for (i=low; i<high; i++) { 3881 if (rp[i] > col) break; 3882 if (rp[i] == col) { 3883 if (is == ADD_VALUES) ap[i] += value; 3884 else ap[i] = value; 3885 goto noinsert; 3886 } 3887 } 3888 if (value == 0.0 && ignorezeroentries) goto noinsert; 3889 if (nonew == 1) goto noinsert; 3890 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3891 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3892 N = nrow++ - 1; a->nz++; high++; 3893 /* shift up all the later entries in this row */ 3894 for (ii=N; ii>=i; ii--) { 3895 rp[ii+1] = rp[ii]; 3896 ap[ii+1] = ap[ii]; 3897 } 3898 rp[i] = col; 3899 ap[i] = value; 3900 noinsert:; 3901 low = i + 1; 3902 } 3903 ailen[row] = nrow; 3904 } 3905 A->same_nonzero = PETSC_FALSE; 3906 PetscFunctionReturnVoid(); 3907 } 3908 EXTERN_C_END 3909