Commits

Hong Zhang committed 0944192

cleanup MatFDColoringApply_SeqBAIJ

  • Participants
  • Parent commits b779938

Comments (0)

Files changed (1)

File src/mat/impls/aij/seq/fdaij.c

 }
                                  
 /* Optimize MatFDColoringApply_AIJ() by using array den2sp to skip calling MatSetValues() */
-/* #define JACOBIANCOLOROPT */
-#if defined(JACOBIANCOLOROPT)
-#include <petsctime.h>
-#endif
 #undef __FUNCT__
 #define __FUNCT__ "MatFDColoringApply_SeqAIJ"
 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
 {
   PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
   PetscErrorCode ierr;
-  PetscInt       k,l,row,col,N,bs=J->rmap->bs,bs2=bs*bs;
+  PetscInt       k,l,row,col,N;
   PetscScalar    dx,*y,*xx,*w3_array;
   PetscScalar    *vscale_array;
   PetscReal      epsilon=coloring->error_rel,umin = coloring->umin,unorm;
   Vec            w1=coloring->w1,w2=coloring->w2,w3;
   void           *fctx=coloring->fctx;
-  PetscBool      flg=PETSC_FALSE,isBAIJ=PETSC_FALSE;
-  PetscScalar    *ca;
+  PetscBool      flg=PETSC_FALSE;
+  Mat_SeqAIJ     *csp=(Mat_SeqAIJ*)J->data;
+  PetscScalar    *ca=csp->a;
   const PetscInt ncolors=coloring->ncolors,*ncolumns=coloring->ncolumns,*nrows=coloring->nrows;
   PetscInt       **columns=coloring->columns,ncolumns_k,nrows_k,nz,spidx;
-#if defined(JACOBIANCOLOROPT)
-  PetscLogDouble t0,t1,t_init=0.0,t_setvals=0.0,t_func=0.0,t_dx=0.0,t_kl=0.0,t00,t11;
-#endif
 
   PetscFunctionBegin;
-#if defined(JACOBIANCOLOROPT)
-  ierr = PetscTime(&t0);CHKERRQ(ierr);
-#endif
-  ierr = PetscObjectTypeCompare((PetscObject)J,MATSEQBAIJ,&isBAIJ);CHKERRQ(ierr);
-  if (isBAIJ) { 
-    Mat_SeqBAIJ     *csp=(Mat_SeqBAIJ*)J->data;
-    ca = csp->a;
-  } else {
-    Mat_SeqAIJ     *csp=(Mat_SeqAIJ*)J->data;
-    ca = csp->a;
-    bs = 1; bs2 = 1;
-  }
-
   ierr = MatSetUnfactored(J);CHKERRQ(ierr);
   ierr = PetscOptionsGetBool(NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);CHKERRQ(ierr);
   if (flg) {
       ierr = MatZeroEntries(J);CHKERRQ(ierr);
     }
   }
-
   if (!coloring->vscale) {
     ierr = VecDuplicate(x1,&coloring->vscale);CHKERRQ(ierr);
   }
 
-  if (coloring->htype[0] == 'w') { /* tacky test; need to make systematic if we add other approaches to computing h*/
-    ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr);
-  }
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t1);CHKERRQ(ierr);
-    t_init += t1 - t0;
-#endif
-
   /* Set w1 = F(x1) */
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t0);CHKERRQ(ierr);
-#endif
   if (!coloring->fset) {
     ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
     ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr);
   } else {
     coloring->fset = PETSC_FALSE;
   }
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t1);CHKERRQ(ierr);
-    t_setvals += t1 - t0;
-#endif
-
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t0);CHKERRQ(ierr);
-#endif
+
   if (!coloring->w3) {
     ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr);
     ierr = PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w3);CHKERRQ(ierr);
   } 
   w3 = coloring->w3;
   
-  /* Compute scale factors: vscale = 1./dx = 1./(epsilon*xx) */
-  ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);
-  ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);
-  ierr = VecGetSize(x1,&N);CHKERRQ(ierr);
-  for (col=0; col<N; col++) {
-    if (coloring->htype[0] == 'w') {
-      dx = 1.0 + unorm;
-    } else {
+  /* Compute scale factors: vscale = dx */
+  if (coloring->htype[0] == 'w') {
+    /* tacky test; need to make systematic if we add other approaches to computing h*/
+    ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr);
+    dx = (1.0 + unorm)*epsilon;
+    ierr = VecSet(coloring->vscale,dx);CHKERRQ(ierr);
+  } else {
+    ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);
+    ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);
+    ierr = VecGetSize(x1,&N);CHKERRQ(ierr);
+    for (col=0; col<N; col++) {
       dx = xx[col];
+      if (dx == 0.0) dx = 1.0;
+      if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
+      else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
+      dx               *= epsilon;
+      vscale_array[col] = dx;
     }
-    if (dx == (PetscScalar)0.0) dx = 1.0;
-    if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0)     dx = umin;
-    else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin;
-    dx               *= epsilon;
-    vscale_array[col] = (PetscScalar)dx;
+    ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);
+    ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr);
   }
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t1);CHKERRQ(ierr);
-    t_init += t1 - t0;
-#endif
  
   nz  = 0;
+  ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);
   for (k=0; k<ncolors; k++) { /* loop over colors */
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t0);CHKERRQ(ierr);
-#endif
     coloring->currentcolor = k;
 
     /*
       w3_array[col] += vscale_array[col]; 
     }
     ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr);
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t1);CHKERRQ(ierr);
-    t_dx += t1 - t0;
-#endif
 
     /*
       Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
                            w2 = F(x1 + dx) - F(x1)
     */
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t0);CHKERRQ(ierr);
-#endif
     ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
     ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr);
     ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t1);CHKERRQ(ierr);
-    t_func += t1 - t0;
-#endif
-
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t0);CHKERRQ(ierr);
-#endif
     ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr);
 
     /*
     ierr = VecGetArray(w2,&y);CHKERRQ(ierr);
     nrows_k = nrows[k];
     for (l=0; l<nrows_k; l++) { /* loop over rows */
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t00);CHKERRQ(ierr);
-#endif
       row   = coloring->rowcolden2sp3[nz++];
       col   = coloring->rowcolden2sp3[nz++];
       spidx = coloring->rowcolden2sp3[nz++];
-#if defined(JACOBIANCOLOROPT)
-      ierr = PetscTime(&t11);CHKERRQ(ierr);
-      t_kl += t11 - t00;
-#endif
-      //printf(" row col val
-      ca[spidx*bs2] = y[row*bs]/vscale_array[col*bs];
-      //printf(" (%d, %d, %g)\n", row*bs,col*bs,ca[spidx*bs2]);
+      ca[spidx] = y[row]/vscale_array[col];
     }
     ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr);
-
-#if defined(JACOBIANCOLOROPT)
-    ierr = PetscTime(&t1);CHKERRQ(ierr);
-    t_setvals += t1 - t0;
-#endif
   } /* endof for each color */
-#if defined(JACOBIANCOLOROPT)
-  printf("     FDColorApply time: init %g + func %g + setvalues %g + dx %g= %g\n",t_init,t_func,t_setvals,t_dx,t_init+t_func+t_setvals+t_dx);
-  printf("     FDColorApply time: kl %g\n",t_kl);
-#endif
   ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);
-  ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr);
 
   coloring->currentcolor = -1;
   PetscFunctionReturn(0);