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