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