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