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