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