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