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