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 = VecGetArrayRead(xx,&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 = VecRestoreArrayRead(xx,&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 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1230 CHKMEMQ; 1231 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1232 if (flag == SOR_APPLY_UPPER) { 1233 /* apply (U + D/omega) to the vector */ 1234 bs = b; 1235 for (i=0; i<m; i++) { 1236 d = fshift + mdiag[i]; 1237 n = a->i[i+1] - diag[i] - 1; 1238 idx = a->j + diag[i] + 1; 1239 v = a->a + diag[i] + 1; 1240 sum = b[i]*d/omega; 1241 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1242 x[i] = sum; 1243 } 1244 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1245 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1246 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1247 PetscFunctionReturn(0); 1248 } 1249 1250 if (flag == SOR_APPLY_LOWER) { 1251 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1252 } else if (flag & SOR_EISENSTAT) { 1253 /* Let A = L + U + D; where L is lower trianglar, 1254 U is upper triangular, E = D/omega; This routine applies 1255 1256 (L + E)^{-1} A (U + E)^{-1} 1257 1258 to a vector efficiently using Eisenstat's trick. 1259 */ 1260 scale = (2.0/omega) - 1.0; 1261 1262 /* x = (E + U)^{-1} b */ 1263 for (i=m-1; i>=0; i--) { 1264 n = a->i[i+1] - diag[i] - 1; 1265 idx = a->j + diag[i] + 1; 1266 v = a->a + diag[i] + 1; 1267 sum = b[i]; 1268 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1269 x[i] = sum*idiag[i]; 1270 } 1271 1272 /* t = b - (2*E - D)x */ 1273 v = a->a; 1274 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1275 1276 /* t = (E + L)^{-1}t */ 1277 ts = t; 1278 diag = a->diag; 1279 for (i=0; i<m; i++) { 1280 n = diag[i] - a->i[i]; 1281 idx = a->j + a->i[i]; 1282 v = a->a + a->i[i]; 1283 sum = t[i]; 1284 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1285 t[i] = sum*idiag[i]; 1286 /* x = x + t */ 1287 x[i] += t[i]; 1288 } 1289 1290 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1291 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1292 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1293 PetscFunctionReturn(0); 1294 } 1295 if (flag & SOR_ZERO_INITIAL_GUESS) { 1296 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1297 for (i=0; i<m; i++) { 1298 n = diag[i] - a->i[i]; 1299 idx = a->j + a->i[i]; 1300 v = a->a + a->i[i]; 1301 sum = b[i]; 1302 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1303 t[i] = sum; 1304 x[i] = sum*idiag[i]; 1305 } 1306 xb = t; 1307 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1308 } else xb = b; 1309 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1310 for (i=m-1; i>=0; i--) { 1311 n = a->i[i+1] - diag[i] - 1; 1312 idx = a->j + diag[i] + 1; 1313 v = a->a + diag[i] + 1; 1314 sum = xb[i]; 1315 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1316 if (xb == b) { 1317 x[i] = sum*idiag[i]; 1318 } else { 1319 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1320 } 1321 } 1322 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1323 } 1324 its--; 1325 } 1326 while (its--) { 1327 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1328 for (i=0; i<m; i++) { 1329 n = a->i[i+1] - a->i[i]; 1330 idx = a->j + a->i[i]; 1331 v = a->a + a->i[i]; 1332 sum = b[i]; 1333 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1334 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1335 } 1336 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1337 } 1338 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1339 for (i=m-1; i>=0; i--) { 1340 n = a->i[i+1] - a->i[i]; 1341 idx = a->j + a->i[i]; 1342 v = a->a + a->i[i]; 1343 sum = b[i]; 1344 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1345 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1346 } 1347 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1348 } 1349 } 1350 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1351 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1352 CHKMEMQ; PetscFunctionReturn(0); 1353 } 1354 1355 1356 #undef __FUNCT__ 1357 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1358 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1359 { 1360 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1361 1362 PetscFunctionBegin; 1363 info->block_size = 1.0; 1364 info->nz_allocated = (double)a->maxnz; 1365 info->nz_used = (double)a->nz; 1366 info->nz_unneeded = (double)(a->maxnz - a->nz); 1367 info->assemblies = (double)A->num_ass; 1368 info->mallocs = (double)A->info.mallocs; 1369 info->memory = ((PetscObject)A)->mem; 1370 if (A->factortype) { 1371 info->fill_ratio_given = A->info.fill_ratio_given; 1372 info->fill_ratio_needed = A->info.fill_ratio_needed; 1373 info->factor_mallocs = A->info.factor_mallocs; 1374 } else { 1375 info->fill_ratio_given = 0; 1376 info->fill_ratio_needed = 0; 1377 info->factor_mallocs = 0; 1378 } 1379 PetscFunctionReturn(0); 1380 } 1381 1382 #undef __FUNCT__ 1383 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1384 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1385 { 1386 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1387 PetscInt i,m = A->rmap->n - 1,d = 0; 1388 PetscErrorCode ierr; 1389 PetscTruth missing; 1390 1391 PetscFunctionBegin; 1392 if (a->keepnonzeropattern) { 1393 for (i=0; i<N; i++) { 1394 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1395 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1396 } 1397 if (diag != 0.0) { 1398 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1399 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1400 for (i=0; i<N; i++) { 1401 a->a[a->diag[rows[i]]] = diag; 1402 } 1403 } 1404 A->same_nonzero = PETSC_TRUE; 1405 } else { 1406 if (diag != 0.0) { 1407 for (i=0; i<N; i++) { 1408 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1409 if (a->ilen[rows[i]] > 0) { 1410 a->ilen[rows[i]] = 1; 1411 a->a[a->i[rows[i]]] = diag; 1412 a->j[a->i[rows[i]]] = rows[i]; 1413 } else { /* in case row was completely empty */ 1414 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1415 } 1416 } 1417 } else { 1418 for (i=0; i<N; i++) { 1419 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1420 a->ilen[rows[i]] = 0; 1421 } 1422 } 1423 A->same_nonzero = PETSC_FALSE; 1424 } 1425 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1426 PetscFunctionReturn(0); 1427 } 1428 1429 #undef __FUNCT__ 1430 #define __FUNCT__ "MatGetRow_SeqAIJ" 1431 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1432 { 1433 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1434 PetscInt *itmp; 1435 1436 PetscFunctionBegin; 1437 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1438 1439 *nz = a->i[row+1] - a->i[row]; 1440 if (v) *v = a->a + a->i[row]; 1441 if (idx) { 1442 itmp = a->j + a->i[row]; 1443 if (*nz) { 1444 *idx = itmp; 1445 } 1446 else *idx = 0; 1447 } 1448 PetscFunctionReturn(0); 1449 } 1450 1451 /* remove this function? */ 1452 #undef __FUNCT__ 1453 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1454 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1455 { 1456 PetscFunctionBegin; 1457 PetscFunctionReturn(0); 1458 } 1459 1460 #undef __FUNCT__ 1461 #define __FUNCT__ "MatNorm_SeqAIJ" 1462 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1463 { 1464 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1465 MatScalar *v = a->a; 1466 PetscReal sum = 0.0; 1467 PetscErrorCode ierr; 1468 PetscInt i,j; 1469 1470 PetscFunctionBegin; 1471 if (type == NORM_FROBENIUS) { 1472 for (i=0; i<a->nz; i++) { 1473 #if defined(PETSC_USE_COMPLEX) 1474 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1475 #else 1476 sum += (*v)*(*v); v++; 1477 #endif 1478 } 1479 *nrm = sqrt(sum); 1480 } else if (type == NORM_1) { 1481 PetscReal *tmp; 1482 PetscInt *jj = a->j; 1483 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1484 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1485 *nrm = 0.0; 1486 for (j=0; j<a->nz; j++) { 1487 tmp[*jj++] += PetscAbsScalar(*v); v++; 1488 } 1489 for (j=0; j<A->cmap->n; j++) { 1490 if (tmp[j] > *nrm) *nrm = tmp[j]; 1491 } 1492 ierr = PetscFree(tmp);CHKERRQ(ierr); 1493 } else if (type == NORM_INFINITY) { 1494 *nrm = 0.0; 1495 for (j=0; j<A->rmap->n; j++) { 1496 v = a->a + a->i[j]; 1497 sum = 0.0; 1498 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1499 sum += PetscAbsScalar(*v); v++; 1500 } 1501 if (sum > *nrm) *nrm = sum; 1502 } 1503 } else { 1504 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1505 } 1506 PetscFunctionReturn(0); 1507 } 1508 1509 #undef __FUNCT__ 1510 #define __FUNCT__ "MatTranspose_SeqAIJ" 1511 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1512 { 1513 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1514 Mat C; 1515 PetscErrorCode ierr; 1516 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1517 MatScalar *array = a->a; 1518 1519 PetscFunctionBegin; 1520 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"); 1521 1522 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1523 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1524 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1525 1526 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1527 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1528 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1529 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1530 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1531 ierr = PetscFree(col);CHKERRQ(ierr); 1532 } else { 1533 C = *B; 1534 } 1535 1536 for (i=0; i<m; i++) { 1537 len = ai[i+1]-ai[i]; 1538 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1539 array += len; 1540 aj += len; 1541 } 1542 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1543 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1544 1545 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1546 *B = C; 1547 } else { 1548 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1549 } 1550 PetscFunctionReturn(0); 1551 } 1552 1553 EXTERN_C_BEGIN 1554 #undef __FUNCT__ 1555 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1556 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1557 { 1558 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1559 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1560 MatScalar *va,*vb; 1561 PetscErrorCode ierr; 1562 PetscInt ma,na,mb,nb, i; 1563 1564 PetscFunctionBegin; 1565 bij = (Mat_SeqAIJ *) B->data; 1566 1567 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1568 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1569 if (ma!=nb || na!=mb){ 1570 *f = PETSC_FALSE; 1571 PetscFunctionReturn(0); 1572 } 1573 aii = aij->i; bii = bij->i; 1574 adx = aij->j; bdx = bij->j; 1575 va = aij->a; vb = bij->a; 1576 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1577 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1578 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1579 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1580 1581 *f = PETSC_TRUE; 1582 for (i=0; i<ma; i++) { 1583 while (aptr[i]<aii[i+1]) { 1584 PetscInt idc,idr; 1585 PetscScalar vc,vr; 1586 /* column/row index/value */ 1587 idc = adx[aptr[i]]; 1588 idr = bdx[bptr[idc]]; 1589 vc = va[aptr[i]]; 1590 vr = vb[bptr[idc]]; 1591 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1592 *f = PETSC_FALSE; 1593 goto done; 1594 } else { 1595 aptr[i]++; 1596 if (B || i!=idc) bptr[idc]++; 1597 } 1598 } 1599 } 1600 done: 1601 ierr = PetscFree(aptr);CHKERRQ(ierr); 1602 if (B) { 1603 ierr = PetscFree(bptr);CHKERRQ(ierr); 1604 } 1605 PetscFunctionReturn(0); 1606 } 1607 EXTERN_C_END 1608 1609 EXTERN_C_BEGIN 1610 #undef __FUNCT__ 1611 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1612 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1613 { 1614 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1615 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1616 MatScalar *va,*vb; 1617 PetscErrorCode ierr; 1618 PetscInt ma,na,mb,nb, i; 1619 1620 PetscFunctionBegin; 1621 bij = (Mat_SeqAIJ *) B->data; 1622 1623 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1624 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1625 if (ma!=nb || na!=mb){ 1626 *f = PETSC_FALSE; 1627 PetscFunctionReturn(0); 1628 } 1629 aii = aij->i; bii = bij->i; 1630 adx = aij->j; bdx = bij->j; 1631 va = aij->a; vb = bij->a; 1632 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1633 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1634 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1635 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1636 1637 *f = PETSC_TRUE; 1638 for (i=0; i<ma; i++) { 1639 while (aptr[i]<aii[i+1]) { 1640 PetscInt idc,idr; 1641 PetscScalar vc,vr; 1642 /* column/row index/value */ 1643 idc = adx[aptr[i]]; 1644 idr = bdx[bptr[idc]]; 1645 vc = va[aptr[i]]; 1646 vr = vb[bptr[idc]]; 1647 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1648 *f = PETSC_FALSE; 1649 goto done; 1650 } else { 1651 aptr[i]++; 1652 if (B || i!=idc) bptr[idc]++; 1653 } 1654 } 1655 } 1656 done: 1657 ierr = PetscFree(aptr);CHKERRQ(ierr); 1658 if (B) { 1659 ierr = PetscFree(bptr);CHKERRQ(ierr); 1660 } 1661 PetscFunctionReturn(0); 1662 } 1663 EXTERN_C_END 1664 1665 #undef __FUNCT__ 1666 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1667 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1668 { 1669 PetscErrorCode ierr; 1670 PetscFunctionBegin; 1671 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1672 PetscFunctionReturn(0); 1673 } 1674 1675 #undef __FUNCT__ 1676 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1677 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1678 { 1679 PetscErrorCode ierr; 1680 PetscFunctionBegin; 1681 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1682 PetscFunctionReturn(0); 1683 } 1684 1685 #undef __FUNCT__ 1686 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1687 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1688 { 1689 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1690 PetscScalar *l,*r,x; 1691 MatScalar *v; 1692 PetscErrorCode ierr; 1693 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1694 1695 PetscFunctionBegin; 1696 if (ll) { 1697 /* The local size is used so that VecMPI can be passed to this routine 1698 by MatDiagonalScale_MPIAIJ */ 1699 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1700 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1701 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1702 v = a->a; 1703 for (i=0; i<m; i++) { 1704 x = l[i]; 1705 M = a->i[i+1] - a->i[i]; 1706 for (j=0; j<M; j++) { (*v++) *= x;} 1707 } 1708 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1709 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1710 } 1711 if (rr) { 1712 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1713 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1714 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1715 v = a->a; jj = a->j; 1716 for (i=0; i<nz; i++) { 1717 (*v++) *= r[*jj++]; 1718 } 1719 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1720 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1721 } 1722 PetscFunctionReturn(0); 1723 } 1724 1725 #undef __FUNCT__ 1726 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1727 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1728 { 1729 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1730 PetscErrorCode ierr; 1731 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1732 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1733 const PetscInt *irow,*icol; 1734 PetscInt nrows,ncols; 1735 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1736 MatScalar *a_new,*mat_a; 1737 Mat C; 1738 PetscTruth stride,sorted; 1739 1740 PetscFunctionBegin; 1741 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1742 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1743 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1744 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1745 1746 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1747 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1748 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1749 1750 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1751 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1752 if (stride && step == 1) { 1753 /* special case of contiguous rows */ 1754 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1755 /* loop over new rows determining lens and starting points */ 1756 for (i=0; i<nrows; i++) { 1757 kstart = ai[irow[i]]; 1758 kend = kstart + ailen[irow[i]]; 1759 for (k=kstart; k<kend; k++) { 1760 if (aj[k] >= first) { 1761 starts[i] = k; 1762 break; 1763 } 1764 } 1765 sum = 0; 1766 while (k < kend) { 1767 if (aj[k++] >= first+ncols) break; 1768 sum++; 1769 } 1770 lens[i] = sum; 1771 } 1772 /* create submatrix */ 1773 if (scall == MAT_REUSE_MATRIX) { 1774 PetscInt n_cols,n_rows; 1775 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1776 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1777 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1778 C = *B; 1779 } else { 1780 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1781 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1782 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1783 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1784 } 1785 c = (Mat_SeqAIJ*)C->data; 1786 1787 /* loop over rows inserting into submatrix */ 1788 a_new = c->a; 1789 j_new = c->j; 1790 i_new = c->i; 1791 1792 for (i=0; i<nrows; i++) { 1793 ii = starts[i]; 1794 lensi = lens[i]; 1795 for (k=0; k<lensi; k++) { 1796 *j_new++ = aj[ii+k] - first; 1797 } 1798 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1799 a_new += lensi; 1800 i_new[i+1] = i_new[i] + lensi; 1801 c->ilen[i] = lensi; 1802 } 1803 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1804 } else { 1805 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1806 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1807 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1808 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1809 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1810 /* determine lens of each row */ 1811 for (i=0; i<nrows; i++) { 1812 kstart = ai[irow[i]]; 1813 kend = kstart + a->ilen[irow[i]]; 1814 lens[i] = 0; 1815 for (k=kstart; k<kend; k++) { 1816 if (smap[aj[k]]) { 1817 lens[i]++; 1818 } 1819 } 1820 } 1821 /* Create and fill new matrix */ 1822 if (scall == MAT_REUSE_MATRIX) { 1823 PetscTruth equal; 1824 1825 c = (Mat_SeqAIJ *)((*B)->data); 1826 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1827 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1828 if (!equal) { 1829 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1830 } 1831 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1832 C = *B; 1833 } else { 1834 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1835 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1836 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1837 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1838 } 1839 c = (Mat_SeqAIJ *)(C->data); 1840 for (i=0; i<nrows; i++) { 1841 row = irow[i]; 1842 kstart = ai[row]; 1843 kend = kstart + a->ilen[row]; 1844 mat_i = c->i[i]; 1845 mat_j = c->j + mat_i; 1846 mat_a = c->a + mat_i; 1847 mat_ilen = c->ilen + i; 1848 for (k=kstart; k<kend; k++) { 1849 if ((tcol=smap[a->j[k]])) { 1850 *mat_j++ = tcol - 1; 1851 *mat_a++ = a->a[k]; 1852 (*mat_ilen)++; 1853 1854 } 1855 } 1856 } 1857 /* Free work space */ 1858 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1859 ierr = PetscFree(smap);CHKERRQ(ierr); 1860 ierr = PetscFree(lens);CHKERRQ(ierr); 1861 } 1862 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1863 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1864 1865 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1866 *B = C; 1867 PetscFunctionReturn(0); 1868 } 1869 1870 #undef __FUNCT__ 1871 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1872 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1873 { 1874 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1875 PetscErrorCode ierr; 1876 Mat outA; 1877 PetscTruth row_identity,col_identity; 1878 1879 PetscFunctionBegin; 1880 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1881 1882 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1883 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1884 1885 outA = inA; 1886 outA->factortype = MAT_FACTOR_LU; 1887 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1888 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1889 a->row = row; 1890 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1891 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1892 a->col = col; 1893 1894 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1895 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1896 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1897 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1898 1899 if (!a->solve_work) { /* this matrix may have been factored before */ 1900 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1901 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1902 } 1903 1904 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1905 if (row_identity && col_identity) { 1906 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1907 } else { 1908 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1909 } 1910 PetscFunctionReturn(0); 1911 } 1912 1913 #undef __FUNCT__ 1914 #define __FUNCT__ "MatScale_SeqAIJ" 1915 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1916 { 1917 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1918 PetscScalar oalpha = alpha; 1919 PetscErrorCode ierr; 1920 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1921 1922 PetscFunctionBegin; 1923 BLASscal_(&bnz,&oalpha,a->a,&one); 1924 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1925 PetscFunctionReturn(0); 1926 } 1927 1928 #undef __FUNCT__ 1929 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1930 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1931 { 1932 PetscErrorCode ierr; 1933 PetscInt i; 1934 1935 PetscFunctionBegin; 1936 if (scall == MAT_INITIAL_MATRIX) { 1937 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1938 } 1939 1940 for (i=0; i<n; i++) { 1941 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1942 } 1943 PetscFunctionReturn(0); 1944 } 1945 1946 #undef __FUNCT__ 1947 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1948 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1949 { 1950 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1951 PetscErrorCode ierr; 1952 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1953 const PetscInt *idx; 1954 PetscInt start,end,*ai,*aj; 1955 PetscBT table; 1956 1957 PetscFunctionBegin; 1958 m = A->rmap->n; 1959 ai = a->i; 1960 aj = a->j; 1961 1962 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1963 1964 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1965 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1966 1967 for (i=0; i<is_max; i++) { 1968 /* Initialize the two local arrays */ 1969 isz = 0; 1970 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1971 1972 /* Extract the indices, assume there can be duplicate entries */ 1973 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1974 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1975 1976 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1977 for (j=0; j<n ; ++j){ 1978 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1979 } 1980 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1981 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1982 1983 k = 0; 1984 for (j=0; j<ov; j++){ /* for each overlap */ 1985 n = isz; 1986 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1987 row = nidx[k]; 1988 start = ai[row]; 1989 end = ai[row+1]; 1990 for (l = start; l<end ; l++){ 1991 val = aj[l] ; 1992 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1993 } 1994 } 1995 } 1996 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1997 } 1998 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1999 ierr = PetscFree(nidx);CHKERRQ(ierr); 2000 PetscFunctionReturn(0); 2001 } 2002 2003 /* -------------------------------------------------------------- */ 2004 #undef __FUNCT__ 2005 #define __FUNCT__ "MatPermute_SeqAIJ" 2006 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2007 { 2008 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2009 PetscErrorCode ierr; 2010 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2011 const PetscInt *row,*col; 2012 PetscInt *cnew,j,*lens; 2013 IS icolp,irowp; 2014 PetscInt *cwork = PETSC_NULL; 2015 PetscScalar *vwork = PETSC_NULL; 2016 2017 PetscFunctionBegin; 2018 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2019 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2020 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2021 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2022 2023 /* determine lengths of permuted rows */ 2024 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2025 for (i=0; i<m; i++) { 2026 lens[row[i]] = a->i[i+1] - a->i[i]; 2027 } 2028 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2029 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2030 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2031 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2032 ierr = PetscFree(lens);CHKERRQ(ierr); 2033 2034 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2035 for (i=0; i<m; i++) { 2036 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2037 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2038 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2039 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2040 } 2041 ierr = PetscFree(cnew);CHKERRQ(ierr); 2042 (*B)->assembled = PETSC_FALSE; 2043 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2044 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2045 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2046 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2047 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2048 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2049 PetscFunctionReturn(0); 2050 } 2051 2052 #undef __FUNCT__ 2053 #define __FUNCT__ "MatCopy_SeqAIJ" 2054 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2055 { 2056 PetscErrorCode ierr; 2057 2058 PetscFunctionBegin; 2059 /* If the two matrices have the same copy implementation, use fast copy. */ 2060 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2061 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2062 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2063 2064 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2065 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2066 } 2067 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2068 } else { 2069 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2070 } 2071 PetscFunctionReturn(0); 2072 } 2073 2074 #undef __FUNCT__ 2075 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2076 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2077 { 2078 PetscErrorCode ierr; 2079 2080 PetscFunctionBegin; 2081 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2082 PetscFunctionReturn(0); 2083 } 2084 2085 #undef __FUNCT__ 2086 #define __FUNCT__ "MatGetArray_SeqAIJ" 2087 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2088 { 2089 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2090 PetscFunctionBegin; 2091 *array = a->a; 2092 PetscFunctionReturn(0); 2093 } 2094 2095 #undef __FUNCT__ 2096 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2097 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2098 { 2099 PetscFunctionBegin; 2100 PetscFunctionReturn(0); 2101 } 2102 2103 #undef __FUNCT__ 2104 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2105 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2106 { 2107 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2108 PetscErrorCode ierr; 2109 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2110 PetscScalar dx,*y,*xx,*w3_array; 2111 PetscScalar *vscale_array; 2112 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2113 Vec w1,w2,w3; 2114 void *fctx = coloring->fctx; 2115 PetscTruth flg = PETSC_FALSE; 2116 2117 PetscFunctionBegin; 2118 if (!coloring->w1) { 2119 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2120 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2121 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2122 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2123 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2124 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2125 } 2126 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2127 2128 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2129 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2130 if (flg) { 2131 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2132 } else { 2133 PetscTruth assembled; 2134 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2135 if (assembled) { 2136 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2137 } 2138 } 2139 2140 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2141 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2142 2143 /* 2144 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2145 coloring->F for the coarser grids from the finest 2146 */ 2147 if (coloring->F) { 2148 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2149 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2150 if (m1 != m2) { 2151 coloring->F = 0; 2152 } 2153 } 2154 2155 if (coloring->F) { 2156 w1 = coloring->F; 2157 coloring->F = 0; 2158 } else { 2159 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2160 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2161 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2162 } 2163 2164 /* 2165 Compute all the scale factors and share with other processors 2166 */ 2167 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2168 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2169 for (k=0; k<coloring->ncolors; k++) { 2170 /* 2171 Loop over each column associated with color adding the 2172 perturbation to the vector w3. 2173 */ 2174 for (l=0; l<coloring->ncolumns[k]; l++) { 2175 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2176 dx = xx[col]; 2177 if (dx == 0.0) dx = 1.0; 2178 #if !defined(PETSC_USE_COMPLEX) 2179 if (dx < umin && dx >= 0.0) dx = umin; 2180 else if (dx < 0.0 && dx > -umin) dx = -umin; 2181 #else 2182 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2183 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2184 #endif 2185 dx *= epsilon; 2186 vscale_array[col] = 1.0/dx; 2187 } 2188 } 2189 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2190 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2191 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2192 2193 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2194 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2195 2196 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2197 else vscaleforrow = coloring->columnsforrow; 2198 2199 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2200 /* 2201 Loop over each color 2202 */ 2203 for (k=0; k<coloring->ncolors; k++) { 2204 coloring->currentcolor = k; 2205 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2206 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2207 /* 2208 Loop over each column associated with color adding the 2209 perturbation to the vector w3. 2210 */ 2211 for (l=0; l<coloring->ncolumns[k]; l++) { 2212 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2213 dx = xx[col]; 2214 if (dx == 0.0) dx = 1.0; 2215 #if !defined(PETSC_USE_COMPLEX) 2216 if (dx < umin && dx >= 0.0) dx = umin; 2217 else if (dx < 0.0 && dx > -umin) dx = -umin; 2218 #else 2219 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2220 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2221 #endif 2222 dx *= epsilon; 2223 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2224 w3_array[col] += dx; 2225 } 2226 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2227 2228 /* 2229 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2230 */ 2231 2232 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2233 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2234 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2235 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2236 2237 /* 2238 Loop over rows of vector, putting results into Jacobian matrix 2239 */ 2240 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2241 for (l=0; l<coloring->nrows[k]; l++) { 2242 row = coloring->rows[k][l]; 2243 col = coloring->columnsforrow[k][l]; 2244 y[row] *= vscale_array[vscaleforrow[k][l]]; 2245 srow = row + start; 2246 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2247 } 2248 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2249 } 2250 coloring->currentcolor = k; 2251 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2252 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2253 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2254 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2255 PetscFunctionReturn(0); 2256 } 2257 2258 #undef __FUNCT__ 2259 #define __FUNCT__ "MatAXPY_SeqAIJ" 2260 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2261 { 2262 PetscErrorCode ierr; 2263 PetscInt i; 2264 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2265 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2266 2267 PetscFunctionBegin; 2268 if (str == SAME_NONZERO_PATTERN) { 2269 PetscScalar alpha = a; 2270 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2271 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2272 if (y->xtoy && y->XtoY != X) { 2273 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2274 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2275 } 2276 if (!y->xtoy) { /* get xtoy */ 2277 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2278 y->XtoY = X; 2279 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2280 } 2281 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2282 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2283 } else { 2284 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2285 } 2286 PetscFunctionReturn(0); 2287 } 2288 2289 #undef __FUNCT__ 2290 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2291 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2292 { 2293 PetscErrorCode ierr; 2294 2295 PetscFunctionBegin; 2296 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2297 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2298 PetscFunctionReturn(0); 2299 } 2300 2301 #undef __FUNCT__ 2302 #define __FUNCT__ "MatConjugate_SeqAIJ" 2303 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2304 { 2305 #if defined(PETSC_USE_COMPLEX) 2306 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2307 PetscInt i,nz; 2308 PetscScalar *a; 2309 2310 PetscFunctionBegin; 2311 nz = aij->nz; 2312 a = aij->a; 2313 for (i=0; i<nz; i++) { 2314 a[i] = PetscConj(a[i]); 2315 } 2316 #else 2317 PetscFunctionBegin; 2318 #endif 2319 PetscFunctionReturn(0); 2320 } 2321 2322 #undef __FUNCT__ 2323 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2324 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2325 { 2326 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2327 PetscErrorCode ierr; 2328 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2329 PetscReal atmp; 2330 PetscScalar *x; 2331 MatScalar *aa; 2332 2333 PetscFunctionBegin; 2334 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2335 aa = a->a; 2336 ai = a->i; 2337 aj = a->j; 2338 2339 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2340 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2341 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2342 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2343 for (i=0; i<m; i++) { 2344 ncols = ai[1] - ai[0]; ai++; 2345 x[i] = 0.0; 2346 for (j=0; j<ncols; j++){ 2347 atmp = PetscAbsScalar(*aa); 2348 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2349 aa++; aj++; 2350 } 2351 } 2352 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2353 PetscFunctionReturn(0); 2354 } 2355 2356 #undef __FUNCT__ 2357 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2358 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2359 { 2360 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2361 PetscErrorCode ierr; 2362 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2363 PetscScalar *x; 2364 MatScalar *aa; 2365 2366 PetscFunctionBegin; 2367 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2368 aa = a->a; 2369 ai = a->i; 2370 aj = a->j; 2371 2372 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2373 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2374 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2375 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2376 for (i=0; i<m; i++) { 2377 ncols = ai[1] - ai[0]; ai++; 2378 if (ncols == A->cmap->n) { /* row is dense */ 2379 x[i] = *aa; if (idx) idx[i] = 0; 2380 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2381 x[i] = 0.0; 2382 if (idx) { 2383 idx[i] = 0; /* in case ncols is zero */ 2384 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2385 if (aj[j] > j) { 2386 idx[i] = j; 2387 break; 2388 } 2389 } 2390 } 2391 } 2392 for (j=0; j<ncols; j++){ 2393 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2394 aa++; aj++; 2395 } 2396 } 2397 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2398 PetscFunctionReturn(0); 2399 } 2400 2401 #undef __FUNCT__ 2402 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2403 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2404 { 2405 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2406 PetscErrorCode ierr; 2407 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2408 PetscReal atmp; 2409 PetscScalar *x; 2410 MatScalar *aa; 2411 2412 PetscFunctionBegin; 2413 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2414 aa = a->a; 2415 ai = a->i; 2416 aj = a->j; 2417 2418 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2419 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2420 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2421 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2422 for (i=0; i<m; i++) { 2423 ncols = ai[1] - ai[0]; ai++; 2424 if (ncols) { 2425 /* Get first nonzero */ 2426 for(j = 0; j < ncols; j++) { 2427 atmp = PetscAbsScalar(aa[j]); 2428 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2429 } 2430 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2431 } else { 2432 x[i] = 0.0; if (idx) idx[i] = 0; 2433 } 2434 for(j = 0; j < ncols; j++) { 2435 atmp = PetscAbsScalar(*aa); 2436 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2437 aa++; aj++; 2438 } 2439 } 2440 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2441 PetscFunctionReturn(0); 2442 } 2443 2444 #undef __FUNCT__ 2445 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2446 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2447 { 2448 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2449 PetscErrorCode ierr; 2450 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2451 PetscScalar *x; 2452 MatScalar *aa; 2453 2454 PetscFunctionBegin; 2455 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2456 aa = a->a; 2457 ai = a->i; 2458 aj = a->j; 2459 2460 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2461 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2462 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2463 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2464 for (i=0; i<m; i++) { 2465 ncols = ai[1] - ai[0]; ai++; 2466 if (ncols == A->cmap->n) { /* row is dense */ 2467 x[i] = *aa; if (idx) idx[i] = 0; 2468 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2469 x[i] = 0.0; 2470 if (idx) { /* find first implicit 0.0 in the row */ 2471 idx[i] = 0; /* in case ncols is zero */ 2472 for (j=0;j<ncols;j++) { 2473 if (aj[j] > j) { 2474 idx[i] = j; 2475 break; 2476 } 2477 } 2478 } 2479 } 2480 for (j=0; j<ncols; j++){ 2481 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2482 aa++; aj++; 2483 } 2484 } 2485 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2486 PetscFunctionReturn(0); 2487 } 2488 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2489 /* -------------------------------------------------------------------*/ 2490 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2491 MatGetRow_SeqAIJ, 2492 MatRestoreRow_SeqAIJ, 2493 MatMult_SeqAIJ, 2494 /* 4*/ MatMultAdd_SeqAIJ, 2495 MatMultTranspose_SeqAIJ, 2496 MatMultTransposeAdd_SeqAIJ, 2497 0, 2498 0, 2499 0, 2500 /*10*/ 0, 2501 MatLUFactor_SeqAIJ, 2502 0, 2503 MatSOR_SeqAIJ, 2504 MatTranspose_SeqAIJ, 2505 /*15*/ MatGetInfo_SeqAIJ, 2506 MatEqual_SeqAIJ, 2507 MatGetDiagonal_SeqAIJ, 2508 MatDiagonalScale_SeqAIJ, 2509 MatNorm_SeqAIJ, 2510 /*20*/ 0, 2511 MatAssemblyEnd_SeqAIJ, 2512 MatSetOption_SeqAIJ, 2513 MatZeroEntries_SeqAIJ, 2514 /*24*/ MatZeroRows_SeqAIJ, 2515 0, 2516 0, 2517 0, 2518 0, 2519 /*29*/ MatSetUpPreallocation_SeqAIJ, 2520 0, 2521 0, 2522 MatGetArray_SeqAIJ, 2523 MatRestoreArray_SeqAIJ, 2524 /*34*/ MatDuplicate_SeqAIJ, 2525 0, 2526 0, 2527 MatILUFactor_SeqAIJ, 2528 0, 2529 /*39*/ MatAXPY_SeqAIJ, 2530 MatGetSubMatrices_SeqAIJ, 2531 MatIncreaseOverlap_SeqAIJ, 2532 MatGetValues_SeqAIJ, 2533 MatCopy_SeqAIJ, 2534 /*44*/ MatGetRowMax_SeqAIJ, 2535 MatScale_SeqAIJ, 2536 0, 2537 MatDiagonalSet_SeqAIJ, 2538 0, 2539 /*49*/ MatSetBlockSize_SeqAIJ, 2540 MatGetRowIJ_SeqAIJ, 2541 MatRestoreRowIJ_SeqAIJ, 2542 MatGetColumnIJ_SeqAIJ, 2543 MatRestoreColumnIJ_SeqAIJ, 2544 /*54*/ MatFDColoringCreate_SeqAIJ, 2545 0, 2546 0, 2547 MatPermute_SeqAIJ, 2548 0, 2549 /*59*/ 0, 2550 MatDestroy_SeqAIJ, 2551 MatView_SeqAIJ, 2552 0, 2553 0, 2554 /*64*/ 0, 2555 0, 2556 0, 2557 0, 2558 0, 2559 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2560 MatGetRowMinAbs_SeqAIJ, 2561 0, 2562 MatSetColoring_SeqAIJ, 2563 #if defined(PETSC_HAVE_ADIC) 2564 MatSetValuesAdic_SeqAIJ, 2565 #else 2566 0, 2567 #endif 2568 /*74*/ MatSetValuesAdifor_SeqAIJ, 2569 MatFDColoringApply_AIJ, 2570 0, 2571 0, 2572 0, 2573 /*79*/ 0, 2574 0, 2575 0, 2576 0, 2577 MatLoad_SeqAIJ, 2578 /*84*/ MatIsSymmetric_SeqAIJ, 2579 MatIsHermitian_SeqAIJ, 2580 0, 2581 0, 2582 0, 2583 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2584 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2585 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2586 MatPtAP_Basic, 2587 MatPtAPSymbolic_SeqAIJ, 2588 /*94*/ MatPtAPNumeric_SeqAIJ, 2589 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2590 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2591 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2592 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2593 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2594 0, 2595 0, 2596 MatConjugate_SeqAIJ, 2597 0, 2598 /*104*/MatSetValuesRow_SeqAIJ, 2599 MatRealPart_SeqAIJ, 2600 MatImaginaryPart_SeqAIJ, 2601 0, 2602 0, 2603 /*109*/0, 2604 0, 2605 MatGetRowMin_SeqAIJ, 2606 0, 2607 MatMissingDiagonal_SeqAIJ, 2608 /*114*/0, 2609 0, 2610 0, 2611 0, 2612 0 2613 }; 2614 2615 EXTERN_C_BEGIN 2616 #undef __FUNCT__ 2617 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2618 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2619 { 2620 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2621 PetscInt i,nz,n; 2622 2623 PetscFunctionBegin; 2624 2625 nz = aij->maxnz; 2626 n = mat->rmap->n; 2627 for (i=0; i<nz; i++) { 2628 aij->j[i] = indices[i]; 2629 } 2630 aij->nz = nz; 2631 for (i=0; i<n; i++) { 2632 aij->ilen[i] = aij->imax[i]; 2633 } 2634 2635 PetscFunctionReturn(0); 2636 } 2637 EXTERN_C_END 2638 2639 #undef __FUNCT__ 2640 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2641 /*@ 2642 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2643 in the matrix. 2644 2645 Input Parameters: 2646 + mat - the SeqAIJ matrix 2647 - indices - the column indices 2648 2649 Level: advanced 2650 2651 Notes: 2652 This can be called if you have precomputed the nonzero structure of the 2653 matrix and want to provide it to the matrix object to improve the performance 2654 of the MatSetValues() operation. 2655 2656 You MUST have set the correct numbers of nonzeros per row in the call to 2657 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2658 2659 MUST be called before any calls to MatSetValues(); 2660 2661 The indices should start with zero, not one. 2662 2663 @*/ 2664 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2665 { 2666 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2667 2668 PetscFunctionBegin; 2669 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2670 PetscValidPointer(indices,2); 2671 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2672 if (f) { 2673 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2674 } else { 2675 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2676 } 2677 PetscFunctionReturn(0); 2678 } 2679 2680 /* ----------------------------------------------------------------------------------------*/ 2681 2682 EXTERN_C_BEGIN 2683 #undef __FUNCT__ 2684 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2685 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2686 { 2687 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2688 PetscErrorCode ierr; 2689 size_t nz = aij->i[mat->rmap->n]; 2690 2691 PetscFunctionBegin; 2692 if (aij->nonew != 1) { 2693 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2694 } 2695 2696 /* allocate space for values if not already there */ 2697 if (!aij->saved_values) { 2698 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2699 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2700 } 2701 2702 /* copy values over */ 2703 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2704 PetscFunctionReturn(0); 2705 } 2706 EXTERN_C_END 2707 2708 #undef __FUNCT__ 2709 #define __FUNCT__ "MatStoreValues" 2710 /*@ 2711 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2712 example, reuse of the linear part of a Jacobian, while recomputing the 2713 nonlinear portion. 2714 2715 Collect on Mat 2716 2717 Input Parameters: 2718 . mat - the matrix (currently only AIJ matrices support this option) 2719 2720 Level: advanced 2721 2722 Common Usage, with SNESSolve(): 2723 $ Create Jacobian matrix 2724 $ Set linear terms into matrix 2725 $ Apply boundary conditions to matrix, at this time matrix must have 2726 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2727 $ boundary conditions again will not change the nonzero structure 2728 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2729 $ ierr = MatStoreValues(mat); 2730 $ Call SNESSetJacobian() with matrix 2731 $ In your Jacobian routine 2732 $ ierr = MatRetrieveValues(mat); 2733 $ Set nonlinear terms in matrix 2734 2735 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2736 $ // build linear portion of Jacobian 2737 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2738 $ ierr = MatStoreValues(mat); 2739 $ loop over nonlinear iterations 2740 $ ierr = MatRetrieveValues(mat); 2741 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2742 $ // call MatAssemblyBegin/End() on matrix 2743 $ Solve linear system with Jacobian 2744 $ endloop 2745 2746 Notes: 2747 Matrix must already be assemblied before calling this routine 2748 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2749 calling this routine. 2750 2751 When this is called multiple times it overwrites the previous set of stored values 2752 and does not allocated additional space. 2753 2754 .seealso: MatRetrieveValues() 2755 2756 @*/ 2757 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2758 { 2759 PetscErrorCode ierr,(*f)(Mat); 2760 2761 PetscFunctionBegin; 2762 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2763 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2764 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2765 2766 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2767 if (f) { 2768 ierr = (*f)(mat);CHKERRQ(ierr); 2769 } else { 2770 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2771 } 2772 PetscFunctionReturn(0); 2773 } 2774 2775 EXTERN_C_BEGIN 2776 #undef __FUNCT__ 2777 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2778 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2779 { 2780 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2781 PetscErrorCode ierr; 2782 PetscInt nz = aij->i[mat->rmap->n]; 2783 2784 PetscFunctionBegin; 2785 if (aij->nonew != 1) { 2786 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2787 } 2788 if (!aij->saved_values) { 2789 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2790 } 2791 /* copy values over */ 2792 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2793 PetscFunctionReturn(0); 2794 } 2795 EXTERN_C_END 2796 2797 #undef __FUNCT__ 2798 #define __FUNCT__ "MatRetrieveValues" 2799 /*@ 2800 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2801 example, reuse of the linear part of a Jacobian, while recomputing the 2802 nonlinear portion. 2803 2804 Collect on Mat 2805 2806 Input Parameters: 2807 . mat - the matrix (currently on AIJ matrices support this option) 2808 2809 Level: advanced 2810 2811 .seealso: MatStoreValues() 2812 2813 @*/ 2814 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2815 { 2816 PetscErrorCode ierr,(*f)(Mat); 2817 2818 PetscFunctionBegin; 2819 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2820 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2821 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2822 2823 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2824 if (f) { 2825 ierr = (*f)(mat);CHKERRQ(ierr); 2826 } else { 2827 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2828 } 2829 PetscFunctionReturn(0); 2830 } 2831 2832 2833 /* --------------------------------------------------------------------------------*/ 2834 #undef __FUNCT__ 2835 #define __FUNCT__ "MatCreateSeqAIJ" 2836 /*@C 2837 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2838 (the default parallel PETSc format). For good matrix assembly performance 2839 the user should preallocate the matrix storage by setting the parameter nz 2840 (or the array nnz). By setting these parameters accurately, performance 2841 during matrix assembly can be increased by more than a factor of 50. 2842 2843 Collective on MPI_Comm 2844 2845 Input Parameters: 2846 + comm - MPI communicator, set to PETSC_COMM_SELF 2847 . m - number of rows 2848 . n - number of columns 2849 . nz - number of nonzeros per row (same for all rows) 2850 - nnz - array containing the number of nonzeros in the various rows 2851 (possibly different for each row) or PETSC_NULL 2852 2853 Output Parameter: 2854 . A - the matrix 2855 2856 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2857 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2858 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2859 2860 Notes: 2861 If nnz is given then nz is ignored 2862 2863 The AIJ format (also called the Yale sparse matrix format or 2864 compressed row storage), is fully compatible with standard Fortran 77 2865 storage. That is, the stored row and column indices can begin at 2866 either one (as in Fortran) or zero. See the users' manual for details. 2867 2868 Specify the preallocated storage with either nz or nnz (not both). 2869 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2870 allocation. For large problems you MUST preallocate memory or you 2871 will get TERRIBLE performance, see the users' manual chapter on matrices. 2872 2873 By default, this format uses inodes (identical nodes) when possible, to 2874 improve numerical efficiency of matrix-vector products and solves. We 2875 search for consecutive rows with the same nonzero structure, thereby 2876 reusing matrix information to achieve increased efficiency. 2877 2878 Options Database Keys: 2879 + -mat_no_inode - Do not use inodes 2880 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2881 - -mat_aij_oneindex - Internally use indexing starting at 1 2882 rather than 0. Note that when calling MatSetValues(), 2883 the user still MUST index entries starting at 0! 2884 2885 Level: intermediate 2886 2887 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2888 2889 @*/ 2890 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2891 { 2892 PetscErrorCode ierr; 2893 2894 PetscFunctionBegin; 2895 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2896 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2897 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2898 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2899 PetscFunctionReturn(0); 2900 } 2901 2902 #undef __FUNCT__ 2903 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2904 /*@C 2905 MatSeqAIJSetPreallocation - For good matrix assembly performance 2906 the user should preallocate the matrix storage by setting the parameter nz 2907 (or the array nnz). By setting these parameters accurately, performance 2908 during matrix assembly can be increased by more than a factor of 50. 2909 2910 Collective on MPI_Comm 2911 2912 Input Parameters: 2913 + B - The matrix-free 2914 . nz - number of nonzeros per row (same for all rows) 2915 - nnz - array containing the number of nonzeros in the various rows 2916 (possibly different for each row) or PETSC_NULL 2917 2918 Notes: 2919 If nnz is given then nz is ignored 2920 2921 The AIJ format (also called the Yale sparse matrix format or 2922 compressed row storage), is fully compatible with standard Fortran 77 2923 storage. That is, the stored row and column indices can begin at 2924 either one (as in Fortran) or zero. See the users' manual for details. 2925 2926 Specify the preallocated storage with either nz or nnz (not both). 2927 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2928 allocation. For large problems you MUST preallocate memory or you 2929 will get TERRIBLE performance, see the users' manual chapter on matrices. 2930 2931 You can call MatGetInfo() to get information on how effective the preallocation was; 2932 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2933 You can also run with the option -info and look for messages with the string 2934 malloc in them to see if additional memory allocation was needed. 2935 2936 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2937 entries or columns indices 2938 2939 By default, this format uses inodes (identical nodes) when possible, to 2940 improve numerical efficiency of matrix-vector products and solves. We 2941 search for consecutive rows with the same nonzero structure, thereby 2942 reusing matrix information to achieve increased efficiency. 2943 2944 Options Database Keys: 2945 + -mat_no_inode - Do not use inodes 2946 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2947 - -mat_aij_oneindex - Internally use indexing starting at 1 2948 rather than 0. Note that when calling MatSetValues(), 2949 the user still MUST index entries starting at 0! 2950 2951 Level: intermediate 2952 2953 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2954 2955 @*/ 2956 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2957 { 2958 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2959 2960 PetscFunctionBegin; 2961 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2962 if (f) { 2963 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2964 } 2965 PetscFunctionReturn(0); 2966 } 2967 2968 EXTERN_C_BEGIN 2969 #undef __FUNCT__ 2970 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2971 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2972 { 2973 Mat_SeqAIJ *b; 2974 PetscTruth skipallocation = PETSC_FALSE; 2975 PetscErrorCode ierr; 2976 PetscInt i; 2977 2978 PetscFunctionBegin; 2979 2980 if (nz == MAT_SKIP_ALLOCATION) { 2981 skipallocation = PETSC_TRUE; 2982 nz = 0; 2983 } 2984 2985 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 2986 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 2987 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2988 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2989 2990 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2991 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2992 if (nnz) { 2993 for (i=0; i<B->rmap->n; i++) { 2994 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]); 2995 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); 2996 } 2997 } 2998 2999 B->preallocated = PETSC_TRUE; 3000 b = (Mat_SeqAIJ*)B->data; 3001 3002 if (!skipallocation) { 3003 if (!b->imax) { 3004 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3005 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3006 } 3007 if (!nnz) { 3008 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3009 else if (nz <= 0) nz = 1; 3010 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3011 nz = nz*B->rmap->n; 3012 } else { 3013 nz = 0; 3014 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3015 } 3016 /* b->ilen will count nonzeros in each row so far. */ 3017 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3018 3019 /* allocate the matrix space */ 3020 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3021 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3022 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3023 b->i[0] = 0; 3024 for (i=1; i<B->rmap->n+1; i++) { 3025 b->i[i] = b->i[i-1] + b->imax[i-1]; 3026 } 3027 b->singlemalloc = PETSC_TRUE; 3028 b->free_a = PETSC_TRUE; 3029 b->free_ij = PETSC_TRUE; 3030 } else { 3031 b->free_a = PETSC_FALSE; 3032 b->free_ij = PETSC_FALSE; 3033 } 3034 3035 b->nz = 0; 3036 b->maxnz = nz; 3037 B->info.nz_unneeded = (double)b->maxnz; 3038 PetscFunctionReturn(0); 3039 } 3040 EXTERN_C_END 3041 3042 #undef __FUNCT__ 3043 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3044 /*@ 3045 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3046 3047 Input Parameters: 3048 + B - the matrix 3049 . i - the indices into j for the start of each row (starts with zero) 3050 . j - the column indices for each row (starts with zero) these must be sorted for each row 3051 - v - optional values in the matrix 3052 3053 Level: developer 3054 3055 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3056 3057 .keywords: matrix, aij, compressed row, sparse, sequential 3058 3059 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3060 @*/ 3061 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3062 { 3063 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3064 PetscErrorCode ierr; 3065 3066 PetscFunctionBegin; 3067 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3068 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3069 if (f) { 3070 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3071 } 3072 PetscFunctionReturn(0); 3073 } 3074 3075 EXTERN_C_BEGIN 3076 #undef __FUNCT__ 3077 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3078 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3079 { 3080 PetscInt i; 3081 PetscInt m,n; 3082 PetscInt nz; 3083 PetscInt *nnz, nz_max = 0; 3084 PetscScalar *values; 3085 PetscErrorCode ierr; 3086 3087 PetscFunctionBegin; 3088 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3089 3090 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3091 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3092 for(i = 0; i < m; i++) { 3093 nz = Ii[i+1]- Ii[i]; 3094 nz_max = PetscMax(nz_max, nz); 3095 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3096 nnz[i] = nz; 3097 } 3098 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3099 ierr = PetscFree(nnz);CHKERRQ(ierr); 3100 3101 if (v) { 3102 values = (PetscScalar*) v; 3103 } else { 3104 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3105 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3106 } 3107 3108 for(i = 0; i < m; i++) { 3109 nz = Ii[i+1] - Ii[i]; 3110 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3111 } 3112 3113 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3114 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3115 3116 if (!v) { 3117 ierr = PetscFree(values);CHKERRQ(ierr); 3118 } 3119 PetscFunctionReturn(0); 3120 } 3121 EXTERN_C_END 3122 3123 #include "../src/mat/impls/dense/seq/dense.h" 3124 #include "private/petscaxpy.h" 3125 3126 #undef __FUNCT__ 3127 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3128 /* 3129 Computes (B'*A')' since computing B*A directly is untenable 3130 3131 n p p 3132 ( ) ( ) ( ) 3133 m ( A ) * n ( B ) = m ( C ) 3134 ( ) ( ) ( ) 3135 3136 */ 3137 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3138 { 3139 PetscErrorCode ierr; 3140 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3141 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3142 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3143 PetscInt i,n,m,q,p; 3144 const PetscInt *ii,*idx; 3145 const PetscScalar *b,*a,*a_q; 3146 PetscScalar *c,*c_q; 3147 3148 PetscFunctionBegin; 3149 m = A->rmap->n; 3150 n = A->cmap->n; 3151 p = B->cmap->n; 3152 a = sub_a->v; 3153 b = sub_b->a; 3154 c = sub_c->v; 3155 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3156 3157 ii = sub_b->i; 3158 idx = sub_b->j; 3159 for (i=0; i<n; i++) { 3160 q = ii[i+1] - ii[i]; 3161 while (q-->0) { 3162 c_q = c + m*(*idx); 3163 a_q = a + m*i; 3164 PetscAXPY(c_q,*b,a_q,m); 3165 idx++; 3166 b++; 3167 } 3168 } 3169 PetscFunctionReturn(0); 3170 } 3171 3172 #undef __FUNCT__ 3173 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3174 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3175 { 3176 PetscErrorCode ierr; 3177 PetscInt m=A->rmap->n,n=B->cmap->n; 3178 Mat Cmat; 3179 3180 PetscFunctionBegin; 3181 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); 3182 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3183 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3184 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3185 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3186 Cmat->assembled = PETSC_TRUE; 3187 *C = Cmat; 3188 PetscFunctionReturn(0); 3189 } 3190 3191 /* ----------------------------------------------------------------*/ 3192 #undef __FUNCT__ 3193 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3194 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3195 { 3196 PetscErrorCode ierr; 3197 3198 PetscFunctionBegin; 3199 if (scall == MAT_INITIAL_MATRIX){ 3200 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3201 } 3202 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3203 PetscFunctionReturn(0); 3204 } 3205 3206 3207 /*MC 3208 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3209 based on compressed sparse row format. 3210 3211 Options Database Keys: 3212 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3213 3214 Level: beginner 3215 3216 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3217 M*/ 3218 3219 EXTERN_C_BEGIN 3220 #if defined(PETSC_HAVE_PASTIX) 3221 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3222 #endif 3223 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3224 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3225 #endif 3226 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3227 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3228 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3229 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3230 #if defined(PETSC_HAVE_MUMPS) 3231 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3232 #endif 3233 #if defined(PETSC_HAVE_SUPERLU) 3234 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3235 #endif 3236 #if defined(PETSC_HAVE_SUPERLU_DIST) 3237 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3238 #endif 3239 #if defined(PETSC_HAVE_SPOOLES) 3240 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3241 #endif 3242 #if defined(PETSC_HAVE_UMFPACK) 3243 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3244 #endif 3245 #if defined(PETSC_HAVE_CHOLMOD) 3246 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3247 #endif 3248 #if defined(PETSC_HAVE_LUSOL) 3249 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3250 #endif 3251 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3252 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3253 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3254 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3255 #endif 3256 EXTERN_C_END 3257 3258 3259 EXTERN_C_BEGIN 3260 #undef __FUNCT__ 3261 #define __FUNCT__ "MatCreate_SeqAIJ" 3262 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3263 { 3264 Mat_SeqAIJ *b; 3265 PetscErrorCode ierr; 3266 PetscMPIInt size; 3267 3268 PetscFunctionBegin; 3269 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3270 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3271 3272 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3273 B->data = (void*)b; 3274 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3275 B->mapping = 0; 3276 b->row = 0; 3277 b->col = 0; 3278 b->icol = 0; 3279 b->reallocs = 0; 3280 b->ignorezeroentries = PETSC_FALSE; 3281 b->roworiented = PETSC_TRUE; 3282 b->nonew = 0; 3283 b->diag = 0; 3284 b->solve_work = 0; 3285 B->spptr = 0; 3286 b->saved_values = 0; 3287 b->idiag = 0; 3288 b->mdiag = 0; 3289 b->ssor_work = 0; 3290 b->omega = 1.0; 3291 b->fshift = 0.0; 3292 b->idiagvalid = PETSC_FALSE; 3293 b->keepnonzeropattern = PETSC_FALSE; 3294 b->xtoy = 0; 3295 b->XtoY = 0; 3296 b->compressedrow.use = PETSC_FALSE; 3297 b->compressedrow.nrows = B->rmap->n; 3298 b->compressedrow.i = PETSC_NULL; 3299 b->compressedrow.rindex = PETSC_NULL; 3300 b->compressedrow.checked = PETSC_FALSE; 3301 B->same_nonzero = PETSC_FALSE; 3302 3303 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3304 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3305 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3306 "MatGetFactor_seqaij_matlab", 3307 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3308 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3309 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3310 #endif 3311 #if defined(PETSC_HAVE_PASTIX) 3312 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3313 "MatGetFactor_seqaij_pastix", 3314 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3315 #endif 3316 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3317 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3318 "MatGetFactor_seqaij_essl", 3319 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3320 #endif 3321 #if defined(PETSC_HAVE_SUPERLU) 3322 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3323 "MatGetFactor_seqaij_superlu", 3324 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3325 #endif 3326 #if defined(PETSC_HAVE_SUPERLU_DIST) 3327 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3328 "MatGetFactor_seqaij_superlu_dist", 3329 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3330 #endif 3331 #if defined(PETSC_HAVE_SPOOLES) 3332 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3333 "MatGetFactor_seqaij_spooles", 3334 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3335 #endif 3336 #if defined(PETSC_HAVE_MUMPS) 3337 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3338 "MatGetFactor_aij_mumps", 3339 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3340 #endif 3341 #if defined(PETSC_HAVE_UMFPACK) 3342 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3343 "MatGetFactor_seqaij_umfpack", 3344 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3345 #endif 3346 #if defined(PETSC_HAVE_CHOLMOD) 3347 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3348 "MatGetFactor_seqaij_cholmod", 3349 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3350 #endif 3351 #if defined(PETSC_HAVE_LUSOL) 3352 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3353 "MatGetFactor_seqaij_lusol", 3354 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3355 #endif 3356 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3357 "MatGetFactor_seqaij_petsc", 3358 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3359 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3360 "MatGetFactorAvailable_seqaij_petsc", 3361 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3362 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3363 "MatGetFactor_seqaij_bas", 3364 MatGetFactor_seqaij_bas); 3365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3366 "MatSeqAIJSetColumnIndices_SeqAIJ", 3367 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3368 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3369 "MatStoreValues_SeqAIJ", 3370 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3371 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3372 "MatRetrieveValues_SeqAIJ", 3373 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3374 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3375 "MatConvert_SeqAIJ_SeqSBAIJ", 3376 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3377 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3378 "MatConvert_SeqAIJ_SeqBAIJ", 3379 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3380 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3381 "MatConvert_SeqAIJ_SeqCSRPERM", 3382 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3383 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3384 "MatConvert_SeqAIJ_SeqCRL", 3385 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3386 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3387 "MatIsTranspose_SeqAIJ", 3388 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3389 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3390 "MatIsHermitianTranspose_SeqAIJ", 3391 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3392 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3393 "MatSeqAIJSetPreallocation_SeqAIJ", 3394 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3395 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3396 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3397 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3398 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3399 "MatReorderForNonzeroDiagonal_SeqAIJ", 3400 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3401 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3402 "MatMatMult_SeqDense_SeqAIJ", 3403 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3404 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3405 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3406 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3407 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3408 "MatMatMultNumeric_SeqDense_SeqAIJ", 3409 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3410 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3411 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3412 PetscFunctionReturn(0); 3413 } 3414 EXTERN_C_END 3415 3416 #undef __FUNCT__ 3417 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3418 /* 3419 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3420 */ 3421 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3422 { 3423 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3424 PetscErrorCode ierr; 3425 PetscInt i,m = A->rmap->n; 3426 3427 PetscFunctionBegin; 3428 c = (Mat_SeqAIJ*)C->data; 3429 3430 C->factortype = A->factortype; 3431 c->row = 0; 3432 c->col = 0; 3433 c->icol = 0; 3434 c->reallocs = 0; 3435 3436 C->assembled = PETSC_TRUE; 3437 3438 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3439 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3440 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3441 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3442 3443 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3444 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3445 for (i=0; i<m; i++) { 3446 c->imax[i] = a->imax[i]; 3447 c->ilen[i] = a->ilen[i]; 3448 } 3449 3450 /* allocate the matrix space */ 3451 if (mallocmatspace){ 3452 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3453 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3454 c->singlemalloc = PETSC_TRUE; 3455 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3456 if (m > 0) { 3457 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3458 if (cpvalues == MAT_COPY_VALUES) { 3459 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3460 } else { 3461 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3462 } 3463 } 3464 } 3465 3466 c->ignorezeroentries = a->ignorezeroentries; 3467 c->roworiented = a->roworiented; 3468 c->nonew = a->nonew; 3469 if (a->diag) { 3470 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3471 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3472 for (i=0; i<m; i++) { 3473 c->diag[i] = a->diag[i]; 3474 } 3475 } else c->diag = 0; 3476 c->solve_work = 0; 3477 c->saved_values = 0; 3478 c->idiag = 0; 3479 c->ssor_work = 0; 3480 c->keepnonzeropattern = a->keepnonzeropattern; 3481 c->free_a = PETSC_TRUE; 3482 c->free_ij = PETSC_TRUE; 3483 c->xtoy = 0; 3484 c->XtoY = 0; 3485 3486 c->nz = a->nz; 3487 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3488 C->preallocated = PETSC_TRUE; 3489 3490 c->compressedrow.use = a->compressedrow.use; 3491 c->compressedrow.nrows = a->compressedrow.nrows; 3492 c->compressedrow.checked = a->compressedrow.checked; 3493 if (a->compressedrow.checked && a->compressedrow.use){ 3494 i = a->compressedrow.nrows; 3495 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3496 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3497 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3498 } else { 3499 c->compressedrow.use = PETSC_FALSE; 3500 c->compressedrow.i = PETSC_NULL; 3501 c->compressedrow.rindex = PETSC_NULL; 3502 } 3503 C->same_nonzero = A->same_nonzero; 3504 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3505 3506 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3507 PetscFunctionReturn(0); 3508 } 3509 3510 #undef __FUNCT__ 3511 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3512 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3513 { 3514 PetscErrorCode ierr; 3515 3516 PetscFunctionBegin; 3517 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3518 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3519 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3520 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3521 PetscFunctionReturn(0); 3522 } 3523 3524 #undef __FUNCT__ 3525 #define __FUNCT__ "MatLoad_SeqAIJ" 3526 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3527 { 3528 Mat_SeqAIJ *a; 3529 Mat B; 3530 PetscErrorCode ierr; 3531 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3532 int fd; 3533 PetscMPIInt size; 3534 MPI_Comm comm; 3535 3536 PetscFunctionBegin; 3537 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3538 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3539 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3540 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3541 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3542 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3543 M = header[1]; N = header[2]; nz = header[3]; 3544 3545 if (nz < 0) { 3546 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3547 } 3548 3549 /* read in row lengths */ 3550 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3551 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3552 3553 /* check if sum of rowlengths is same as nz */ 3554 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3555 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); 3556 3557 /* create our matrix */ 3558 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3559 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3560 ierr = MatSetType(B,type);CHKERRQ(ierr); 3561 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3562 a = (Mat_SeqAIJ*)B->data; 3563 3564 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3565 3566 /* read in nonzero values */ 3567 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3568 3569 /* set matrix "i" values */ 3570 a->i[0] = 0; 3571 for (i=1; i<= M; i++) { 3572 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3573 a->ilen[i-1] = rowlengths[i-1]; 3574 } 3575 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3576 3577 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3578 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3579 *A = B; 3580 PetscFunctionReturn(0); 3581 } 3582 3583 #undef __FUNCT__ 3584 #define __FUNCT__ "MatEqual_SeqAIJ" 3585 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3586 { 3587 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3588 PetscErrorCode ierr; 3589 #if defined(PETSC_USE_COMPLEX) 3590 PetscInt k; 3591 #endif 3592 3593 PetscFunctionBegin; 3594 /* If the matrix dimensions are not equal,or no of nonzeros */ 3595 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3596 *flg = PETSC_FALSE; 3597 PetscFunctionReturn(0); 3598 } 3599 3600 /* if the a->i are the same */ 3601 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3602 if (!*flg) PetscFunctionReturn(0); 3603 3604 /* if a->j are the same */ 3605 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3606 if (!*flg) PetscFunctionReturn(0); 3607 3608 /* if a->a are the same */ 3609 #if defined(PETSC_USE_COMPLEX) 3610 for (k=0; k<a->nz; k++){ 3611 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3612 *flg = PETSC_FALSE; 3613 PetscFunctionReturn(0); 3614 } 3615 } 3616 #else 3617 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3618 #endif 3619 PetscFunctionReturn(0); 3620 } 3621 3622 #undef __FUNCT__ 3623 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3624 /*@ 3625 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3626 provided by the user. 3627 3628 Collective on MPI_Comm 3629 3630 Input Parameters: 3631 + comm - must be an MPI communicator of size 1 3632 . m - number of rows 3633 . n - number of columns 3634 . i - row indices 3635 . j - column indices 3636 - a - matrix values 3637 3638 Output Parameter: 3639 . mat - the matrix 3640 3641 Level: intermediate 3642 3643 Notes: 3644 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3645 once the matrix is destroyed 3646 3647 You cannot set new nonzero locations into this matrix, that will generate an error. 3648 3649 The i and j indices are 0 based 3650 3651 The format which is used for the sparse matrix input, is equivalent to a 3652 row-major ordering.. i.e for the following matrix, the input data expected is 3653 as shown: 3654 3655 1 0 0 3656 2 0 3 3657 4 5 6 3658 3659 i = {0,1,3,6} [size = nrow+1 = 3+1] 3660 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3661 v = {1,2,3,4,5,6} [size = nz = 6] 3662 3663 3664 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3665 3666 @*/ 3667 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3668 { 3669 PetscErrorCode ierr; 3670 PetscInt ii; 3671 Mat_SeqAIJ *aij; 3672 #if defined(PETSC_USE_DEBUG) 3673 PetscInt jj; 3674 #endif 3675 3676 PetscFunctionBegin; 3677 if (i[0]) { 3678 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3679 } 3680 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3681 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3682 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3683 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3684 aij = (Mat_SeqAIJ*)(*mat)->data; 3685 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3686 3687 aij->i = i; 3688 aij->j = j; 3689 aij->a = a; 3690 aij->singlemalloc = PETSC_FALSE; 3691 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3692 aij->free_a = PETSC_FALSE; 3693 aij->free_ij = PETSC_FALSE; 3694 3695 for (ii=0; ii<m; ii++) { 3696 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3697 #if defined(PETSC_USE_DEBUG) 3698 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]); 3699 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3700 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); 3701 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); 3702 } 3703 #endif 3704 } 3705 #if defined(PETSC_USE_DEBUG) 3706 for (ii=0; ii<aij->i[m]; ii++) { 3707 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3708 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]); 3709 } 3710 #endif 3711 3712 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3713 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3714 PetscFunctionReturn(0); 3715 } 3716 3717 #undef __FUNCT__ 3718 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3719 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3720 { 3721 PetscErrorCode ierr; 3722 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3723 3724 PetscFunctionBegin; 3725 if (coloring->ctype == IS_COLORING_GLOBAL) { 3726 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3727 a->coloring = coloring; 3728 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3729 PetscInt i,*larray; 3730 ISColoring ocoloring; 3731 ISColoringValue *colors; 3732 3733 /* set coloring for diagonal portion */ 3734 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3735 for (i=0; i<A->cmap->n; i++) { 3736 larray[i] = i; 3737 } 3738 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3739 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3740 for (i=0; i<A->cmap->n; i++) { 3741 colors[i] = coloring->colors[larray[i]]; 3742 } 3743 ierr = PetscFree(larray);CHKERRQ(ierr); 3744 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3745 a->coloring = ocoloring; 3746 } 3747 PetscFunctionReturn(0); 3748 } 3749 3750 #if defined(PETSC_HAVE_ADIC) 3751 EXTERN_C_BEGIN 3752 #include "adic/ad_utils.h" 3753 EXTERN_C_END 3754 3755 #undef __FUNCT__ 3756 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3757 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3758 { 3759 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3760 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3761 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3762 ISColoringValue *color; 3763 3764 PetscFunctionBegin; 3765 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3766 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3767 color = a->coloring->colors; 3768 /* loop over rows */ 3769 for (i=0; i<m; i++) { 3770 nz = ii[i+1] - ii[i]; 3771 /* loop over columns putting computed value into matrix */ 3772 for (j=0; j<nz; j++) { 3773 *v++ = values[color[*jj++]]; 3774 } 3775 values += nlen; /* jump to next row of derivatives */ 3776 } 3777 PetscFunctionReturn(0); 3778 } 3779 #endif 3780 3781 #undef __FUNCT__ 3782 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3783 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3784 { 3785 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3786 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3787 MatScalar *v = a->a; 3788 PetscScalar *values = (PetscScalar *)advalues; 3789 ISColoringValue *color; 3790 3791 PetscFunctionBegin; 3792 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3793 color = a->coloring->colors; 3794 /* loop over rows */ 3795 for (i=0; i<m; i++) { 3796 nz = ii[i+1] - ii[i]; 3797 /* loop over columns putting computed value into matrix */ 3798 for (j=0; j<nz; j++) { 3799 *v++ = values[color[*jj++]]; 3800 } 3801 values += nl; /* jump to next row of derivatives */ 3802 } 3803 PetscFunctionReturn(0); 3804 } 3805 3806 /* 3807 Special version for direct calls from Fortran 3808 */ 3809 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3810 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3811 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3812 #define matsetvaluesseqaij_ matsetvaluesseqaij 3813 #endif 3814 3815 /* Change these macros so can be used in void function */ 3816 #undef CHKERRQ 3817 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3818 #undef SETERRQ2 3819 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3820 3821 EXTERN_C_BEGIN 3822 #undef __FUNCT__ 3823 #define __FUNCT__ "matsetvaluesseqaij_" 3824 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3825 { 3826 Mat A = *AA; 3827 PetscInt m = *mm, n = *nn; 3828 InsertMode is = *isis; 3829 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3830 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3831 PetscInt *imax,*ai,*ailen; 3832 PetscErrorCode ierr; 3833 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3834 MatScalar *ap,value,*aa; 3835 PetscTruth ignorezeroentries = a->ignorezeroentries; 3836 PetscTruth roworiented = a->roworiented; 3837 3838 PetscFunctionBegin; 3839 ierr = MatPreallocated(A);CHKERRQ(ierr); 3840 imax = a->imax; 3841 ai = a->i; 3842 ailen = a->ilen; 3843 aj = a->j; 3844 aa = a->a; 3845 3846 for (k=0; k<m; k++) { /* loop over added rows */ 3847 row = im[k]; 3848 if (row < 0) continue; 3849 #if defined(PETSC_USE_DEBUG) 3850 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3851 #endif 3852 rp = aj + ai[row]; ap = aa + ai[row]; 3853 rmax = imax[row]; nrow = ailen[row]; 3854 low = 0; 3855 high = nrow; 3856 for (l=0; l<n; l++) { /* loop over added columns */ 3857 if (in[l] < 0) continue; 3858 #if defined(PETSC_USE_DEBUG) 3859 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3860 #endif 3861 col = in[l]; 3862 if (roworiented) { 3863 value = v[l + k*n]; 3864 } else { 3865 value = v[k + l*m]; 3866 } 3867 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3868 3869 if (col <= lastcol) low = 0; else high = nrow; 3870 lastcol = col; 3871 while (high-low > 5) { 3872 t = (low+high)/2; 3873 if (rp[t] > col) high = t; 3874 else low = t; 3875 } 3876 for (i=low; i<high; i++) { 3877 if (rp[i] > col) break; 3878 if (rp[i] == col) { 3879 if (is == ADD_VALUES) ap[i] += value; 3880 else ap[i] = value; 3881 goto noinsert; 3882 } 3883 } 3884 if (value == 0.0 && ignorezeroentries) goto noinsert; 3885 if (nonew == 1) goto noinsert; 3886 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3887 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3888 N = nrow++ - 1; a->nz++; high++; 3889 /* shift up all the later entries in this row */ 3890 for (ii=N; ii>=i; ii--) { 3891 rp[ii+1] = rp[ii]; 3892 ap[ii+1] = ap[ii]; 3893 } 3894 rp[i] = col; 3895 ap[i] = value; 3896 noinsert:; 3897 low = i + 1; 3898 } 3899 ailen[row] = nrow; 3900 } 3901 A->same_nonzero = PETSC_FALSE; 3902 PetscFunctionReturnVoid(); 3903 } 3904 EXTERN_C_END 3905