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