Add ability to apply left- and right-scaling to the sparse Jacobian.

The scaling was optimized in the 1D case to minimize the condition
number of the preconditioned spectral operator.  Unfortunately, it does
not seem to be a clear win on real problems (it is usually worse).  It
can be activated by -djac_tensor_{mass,laplace}_scale TRUE if supported
by the assembly function (see ellip.c for example).

include/dohpfs.h

@@ -137,6 +137,16 @@ EXTERN dErr dFSGetGeometryVector(dFS,Vec*);
                                          {(nx)[(c)[6]][0],(nx)[(c)[6]][1],(nx)[(c)[6]][2]}, \
                                          {(nx)[(c)[7]][0],(nx)[(c)[7]][1],(nx)[(c)[7]][2]}}
 
+#define dQ1SCALE_DECLARE(tscale,scale,i,j,k)                            \
+  const dReal dUNUSED (scale)[8] = {(tscale)[0][(i)+0]*(tscale)[1][(j)+0]*(tscale)[2][(k)+0], \
+                                    (tscale)[0][(i)+1]*(tscale)[1][(j)+0]*(tscale)[2][(k)+0], \
+                                    (tscale)[0][(i)+1]*(tscale)[1][(j)+1]*(tscale)[2][(k)+0], \
+                                    (tscale)[0][(i)+0]*(tscale)[1][(j)+1]*(tscale)[2][(k)+0], \
+                                    (tscale)[0][(i)+0]*(tscale)[1][(j)+0]*(tscale)[2][(k)+1], \
+                                    (tscale)[0][(i)+1]*(tscale)[1][(j)+0]*(tscale)[2][(k)+1], \
+                                    (tscale)[0][(i)+1]*(tscale)[1][(j)+1]*(tscale)[2][(k)+1], \
+                                    (tscale)[0][(i)+0]*(tscale)[1][(j)+1]*(tscale)[2][(k)+1]}
+
 EXTERN dErr dQ1HexComputeQuadrature(const dReal x[8][3],dInt *n,const dReal (**qx)[3],const dReal **jw,const dReal **basis,const dReal **deriv);
 
 

include/dohpjacobi.h

@@ -133,7 +133,7 @@ EXTERN dErr dRuleComputeGeometry(dRule rule,const dReal vtx[restrict][3],dReal[r
 
 EXTERN dErr dEFSView(dEFS efs,PetscViewer viewer);
 EXTERN dErr dEFSGetSizes(dEFS efs,dInt*,dInt *inodes,dInt *total);
-EXTERN dErr dEFSGetTensorNodes(dEFS,dInt*,dInt*,dReal**);
+EXTERN dErr dEFSGetTensorNodes(dEFS,dInt*,dInt*,dReal**,const dReal**,const dReal**);
 EXTERN dErr dEFSGetGlobalCoordinates(dEFS,const dReal vtx[restrict][3],dInt*,dInt[3],dReal(*)[3]);
 EXTERN dErr dEFSGetRule(dEFS efs,dRule *rule);
 EXTERN dErr dEFSApply(dEFS,const dReal[],dInt,const dScalar[],dScalar[restrict],dApplyMode,InsertMode);

include/private/jacimpl.h

@@ -29,7 +29,7 @@ struct _dRuleOps {
 struct _dEFSOps {
   dErr (*view)(dEFS,PetscViewer);
   dErr (*getSizes)(dEFS,dInt*,dInt*,dInt*); /**< topological dimension, number of interior nodes, total number of nodes */
-  dErr (*getTensorNodes)(dEFS,dInt*,dInt*,dReal**);
+  dErr (*getTensorNodes)(dEFS,dInt*,dInt*,dReal**,const dReal**,const dReal**);
   dErr (*apply)(dEFS,const dReal[],dInt,const dScalar[],dScalar[],dApplyMode,InsertMode);
   dErr (*getGlobalCoordinates)(dEFS,const dReal(*)[3],dInt*,dInt[],dReal(*)[3]);
   /**< dofs/node, work length, work, modal values, nodal values */

src/fs/mesh/interface/mesh.c

@@ -534,7 +534,8 @@ dErr dMeshGetEntsOff(dMesh mesh,dMeshESH set,dInt toff[],dMeshEH **inents)
   dFunctionBegin;
   dValidHeader(mesh,dMESH_COOKIE,1);
   dValidPointer(toff,3);
-  dValidPointer(ents,4);
+  dValidPointer(inents,4);
+  *inents = 0;
   err = dMeshGetNumEnts(mesh,set,dTYPE_ALL,dTOPO_ALL,&n);dCHK(err);
   err = dMallocA(n,&ents);dCHK(err);
   toff[0] = 0;

src/fs/tests/ellip.c

@@ -9,6 +9,7 @@ static const char help[] = "Solve a scalar elliptic problem, a regularized p-Lap
 #include "dohpvec.h"
 #include "petscsnes.h"
 
+PetscLogEvent LOG_EllipShellMatMult;
 
 struct EllipParam {
   dReal epsilon;
@@ -352,6 +353,7 @@ static dErr EllipShellMatMult(Mat J,Vec gx,Vec gy)
   dErr err;
 
   dFunctionBegin;
+  err = PetscLogEventBegin(LOG_EllipShellMatMult,J,0,0,0);dCHK(err);
   err = MatShellGetContext(J,(void**)&elp);dCHK(err);
   fs = elp->fs;
   err = dFSGlobalToExpanded(fs,gx,elp->x,dFS_HOMOGENEOUS,INSERT_VALUES);dCHK(err);
@@ -378,6 +380,7 @@ static dErr EllipShellMatMult(Mat J,Vec gx,Vec gy)
   err = VecRestoreArray(elp->x,&x);dCHK(err);
   err = VecRestoreArray(elp->y,&y);dCHK(err);
   err = dFSExpandedToGlobal(fs,elp->y,gy,dFS_HOMOGENEOUS,ADD_VALUES);dCHK(err);
+  err = PetscLogEventEnd(LOG_EllipShellMatMult,J,0,0,0);dCHK(err);
   dFunctionReturn(0);
 }
 
@@ -401,12 +404,16 @@ static dErr EllipJacobian(SNES dUNUSED snes,Vec gx,Mat *J,Mat *Jp,MatStructure *
   err = dFSGetWorkspace(fs,__func__,&nx,NULL,NULL,NULL,NULL,NULL,NULL);dCHK(err); /* We only need space for nodal coordinates */
   for (dInt e=0; e<n; e++) {
     dInt three,P[3];
+    const dReal *tmscale[3],*tlscale[3];
     err = dEFSGetGlobalCoordinates(&efs[e],(const dReal(*)[3])(geom+geomoff[e]),&three,P,nx);dCHK(err);
+    err = dEFSGetTensorNodes(&efs[e],NULL,NULL,NULL,tmscale,tlscale);dCHK(err);
     if (three != 3) dERROR(1,"Dimension not equal to 3");
     for (dInt i=0; i<P[0]-1; i++) { /* P-1 = number of sub-elements in each direction */
       for (dInt j=0; j<P[1]-1; j++) {
         for (dInt k=0; k<P[2]-1; k++) {
           dQ1CORNER_CONST_DECLARE(c,rowcol,corners,off[e],nx,P,i,j,k);
+          dQ1SCALE_DECLARE(tmscale,mscale,i,j,k);
+          dQ1SCALE_DECLARE(tlscale,lscale,i,j,k);
           const dScalar (*uc)[1] = (const dScalar(*)[1])x+off[e]; /* function values, indexed at subelement corners \c uc[c[#]][0] */
           const dReal (*qx)[3],*jw,(*basis)[8],(*deriv)[8][3];
           dInt qn;
@@ -430,14 +437,10 @@ static dErr EllipJacobian(SNES dUNUSED snes,Vec gx,Mat *J,Mat *Jp,MatStructure *
                 const dReal *u = &basis[lq][lp],*Du = deriv[lq][lp];
                 dReal v[1],Dv[3];
                 EllipPointwiseJacobian(&elp->param,&st,jw[lq],u,Du,v,Dv);
-#if 1
                 K[ltest][lp] += //basis[lq][ltest] * v[0]
-                  + deriv[lq][ltest][0] * Dv[0]
-                  + deriv[lq][ltest][1] * Dv[1]
-                  + deriv[lq][ltest][2] * Dv[2];
-#else
-                K[ltest][lp] += basis[lq][ltest] * jw[lq] * basis[lq][lp];
-#endif
+                  lscale[ltest]*lscale[lp]*(+ deriv[lq][ltest][0] * Dv[0]
+                                            + deriv[lq][ltest][1] * Dv[1]
+                                            + deriv[lq][ltest][2] * Dv[2]);
               }
             }
           }
@@ -454,7 +457,7 @@ static dErr EllipJacobian(SNES dUNUSED snes,Vec gx,Mat *J,Mat *Jp,MatStructure *
   err = MatAssemblyEnd(*Jp,MAT_FINAL_ASSEMBLY);dCHK(err);
   err = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);dCHK(err);
   err = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);dCHK(err);
-  *structure = DIFFERENT_NONZERO_PATTERN;
+  *structure = SAME_NONZERO_PATTERN;
   dFunctionReturn(0);
 }
 
@@ -562,6 +565,8 @@ int main(int argc,char *argv[])
   comm = PETSC_COMM_WORLD;
   viewer = PETSC_VIEWER_STDOUT_WORLD;
 
+  err = PetscLogEventRegister("EllipShellMult",MAT_COOKIE,&LOG_EllipShellMatMult);dCHK(err);
+
   err = EllipCreate(comm,&elp);dCHK(err);
   err = EllipSetFromOptions(elp);dCHK(err);
   fs = elp->fs;

src/jacobi/impls/tensor/efstopo.c

@@ -16,7 +16,7 @@ _F(dEFSGetSizes_Tensor_Line);
 _F(dEFSGetSizes_Tensor_Quad);
 _F(dEFSGetSizes_Tensor_Hex);
 #undef _F
-#define _F(f) static dErr f(dEFS,dInt*,dInt*,dReal**restrict) /* dEFSGetTensorNodes */
+#define _F(f) static dErr f(dEFS,dInt*,dInt*,dReal**,const dReal**,const dReal**) /* dEFSGetTensorNodes */
 _F(dEFSGetTensorNodes_Tensor_Line);
 _F(dEFSGetTensorNodes_Tensor_Quad);
 _F(dEFSGetTensorNodes_Tensor_Hex);
@@ -195,7 +195,7 @@ static dErr dEFSGetSizes_Tensor_Hex(dEFS efs,dInt *dim,dInt *inodes,dInt *total)
   dFunctionReturn(0);
 }
 
-static dErr dEFSGetTensorNodes_Tensor_Line(dEFS efs,dInt *dim,dInt tsize[restrict],dReal *x[restrict])
+static dErr dEFSGetTensorNodes_Tensor_Line(dEFS efs,dInt *dim,dInt tsize[],dReal *x[],const dReal *mscale[],const dReal *lscale[])
 {
   TensorBasis *b = ((dEFS_Tensor*)efs)->basis;
 
@@ -211,10 +211,20 @@ static dErr dEFSGetTensorNodes_Tensor_Line(dEFS efs,dInt *dim,dInt tsize[restric
     x[1] = NULL;
     x[2] = NULL;
   }
+  if (mscale) {
+    mscale[0] = b[0]->mscale;
+    mscale[1] = NULL;
+    mscale[2] = NULL;
+  }
+  if (lscale) {
+    lscale[0] = b[0]->lscale;
+    lscale[1] = NULL;
+    lscale[2] = NULL;
+  }
   dFunctionReturn(0);
 }
 
-static dErr dEFSGetTensorNodes_Tensor_Quad(dEFS efs,dInt *dim,dInt tsize[restrict],dReal *x[restrict])
+static dErr dEFSGetTensorNodes_Tensor_Quad(dEFS efs,dInt *dim,dInt tsize[],dReal *x[],const dReal *mscale[],const dReal *lscale[])
 {
   TensorBasis *b = ((dEFS_Tensor*)efs)->basis;
 
@@ -230,10 +240,20 @@ static dErr dEFSGetTensorNodes_Tensor_Quad(dEFS efs,dInt *dim,dInt tsize[restric
     x[1] = b[1]->node;
     x[2] = NULL;
   }
+  if (mscale) {
+    mscale[0] = b[0]->mscale;
+    mscale[1] = b[1]->mscale;
+    mscale[2] = NULL;
+  }
+  if (lscale) {
+    lscale[0] = b[0]->lscale;
+    lscale[1] = b[1]->lscale;
+    lscale[2] = NULL;
+  }
   dFunctionReturn(0);
 }
 
-static dErr dEFSGetTensorNodes_Tensor_Hex(dEFS efs,dInt *dim,dInt tsize[restrict],dReal *x[restrict])
+static dErr dEFSGetTensorNodes_Tensor_Hex(dEFS efs,dInt *dim,dInt tsize[],dReal *x[],const dReal *mscale[],const dReal *lscale[])
 {
   TensorBasis *b = ((dEFS_Tensor*)efs)->basis;
 
@@ -249,6 +269,16 @@ static dErr dEFSGetTensorNodes_Tensor_Hex(dEFS efs,dInt *dim,dInt tsize[restrict
     x[1] = b[1]->node;
     x[2] = b[2]->node;
   }
+  if (mscale) {
+    mscale[0] = b[0]->mscale;
+    mscale[1] = b[1]->mscale;
+    mscale[2] = b[2]->lscale;
+  }
+  if (lscale) {
+    lscale[0] = b[0]->lscale;
+    lscale[1] = b[1]->lscale;
+    lscale[2] = b[2]->lscale;
+  }
   dFunctionReturn(0);
 }
 

src/jacobi/impls/tensor/optimalscale.h

@@ -0,0 +1,53 @@
+/**
+* @file   optimalscale.h
+* @author Jed Brown <jed@59A2.org>
+* @date   Sun May 10 14:18:01 2009
+* 
+* @brief  Defines optimal scaling for q1 matrices assembled on Legendre-Gauss_Lobatto points.
+*
+* These arrays were generated by optimizing over all scalings of the form
+*
+* As = diag(s) * A * diag(s)
+*
+* where A is either the Q1 mass or Laplacian matrix, to minimize the "effective condition number" of inv(As)*A.  The
+* "effective condition number" is computed from the SVD by dropping singular values that are 0 (and would be supplied by
+* boundary values).
+*
+**/
+
+#if !defined(_OPTIMALSCALE_H)
+#define _OPTIMALSCALE_H
+
+#include "dohptype.h"
+
+static const dReal optimal_ones[] = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
+
+static const dReal optimal_mscale_2[] = {1,1},
+  optimal_mscale_3[] = {0.8633789062499998, 1.2732421875000004, 0.8633789062499998},
+  optimal_mscale_4[] = {0.87558593749999991, 1.1244140625000001, 1.1244140625000001, 0.87558593749999991},
+  optimal_mscale_5[] = {0.84529011510312568, 1.1746890407055619, 0.96004168838262516, 1.1746890407055619, 0.84529011510312568},
+  optimal_mscale_6[] = {0.83414557335072448, 1.1347409904428778, 1.0311134362063976, 1.0311134362063976, 1.1347409904428778, 0.83414557335072448},
+  optimal_mscale_7[] = {0.78555955394946975, 1.1276363257065061, 1.0360325284064045, 1.1015431838752399, 1.0360325284064045, 1.1276363257065061, 0.78555955394946975},
+  optimal_mscale_8[] = {0.78005408267437315, 1.115947360985226, 1.0720759672236908, 1.0319225891167099, 1.0319225891167099, 1.0720759672236908, 1.115947360985226, 0.78005408267437315},
+  optimal_mscale_9[] = {0.76452030781575819, 1.0943844400157769, 1.065753105798533, 1.1126528480337201, 0.92537859667242406, 1.1126528480337201, 1.065753105798533, 1.0943844400157769, 0.76452030781575819},
+  optimal_mscale_10[] = {0.76899262661096279, 1.0844808331846316, 1.0650732246613432, 1.1092630792818661, 0.97219023626119672, 0.97219023626119672, 1.1092630792818661, 1.0650732246613432, 1.0844808331846316, 0.76899262661096279},
+  optimal_mscale_11[] = {0.7359699636059529, 1.0641021099995149, 1.0486739809231413, 1.0951863533499089, 1.0124160181349324, 1.0873031479730995, 1.0124160181349324, 1.0951863533499089, 1.0486739809231413, 1.0641021099995149, 0.7359699636059529},
+  optimal_mscale_12[] = {0.75147025229575348, 1.0656982122596839, 1.0354804945478802, 1.0969247634799855, 1.0244541215961656, 1.0259721558205308, 1.0259721558205308, 1.0244541215961656, 1.0969247634799855, 1.0354804945478802, 1.0656982122596839, 0.75147025229575348},
+  optimal_mscale_13[] = {0.73513052753657493, 1.0549675638391163, 1.0187361597988762, 1.0881864667657033, 1.0221852293656799, 1.0598735701897892, 1.0418409650085196, 1.0598735701897892, 1.0221852293656799, 1.0881864667657033, 1.0187361597988762, 1.0549675638391163, 0.73513052753657493},
+  optimal_mscale_14[] = {0.74324111138944104, 1.0579035612110257, 0.99754899979983325, 1.0843431193507969, 1.0105256581725728, 1.0340452975350765, 1.0723922525412544, 1.0723922525412544, 1.0340452975350765, 1.0105256581725728, 1.0843431193507969, 0.99754899979983325, 1.0579035612110257, 0.74324111138944104};
+
+static const dReal optimal_lscale_2[] = {1,1},
+  optimal_lscale_3[] = {0.93129882812499987, 1.1374023437500003, 0.93129882812499987},
+  optimal_lscale_4[] = {0.95996093749999978, 1.0400390625000002, 1.0400390625000002, 0.95996093749999978},
+  optimal_lscale_5[] = {0.95323834643054373, 1.0013589126035178, 1.0908054819318769, 1.0013589126035178, 0.95323834643054373},
+  optimal_lscale_6[] = {0.97233762780379052, 0.99478205400373554, 1.0328803181924742, 1.0328803181924742, 0.99478205400373554, 0.97233762780379052},
+  optimal_lscale_7[] = {0.97491271419538439, 0.9875817846736441, 1.0112409122406048, 1.0525291777807337, 1.0112409122406048, 0.9875817846736441, 0.97491271419538439},
+  optimal_lscale_8[] = {0.98242904574845613, 0.98945640534627399, 1.0008163977013877, 1.0272981512038823, 1.0272981512038823, 1.0008163977013877, 0.98945640534627399, 0.98242904574845613},
+  optimal_lscale_9[] = {0.98371234834817733, 0.98934385044409012, 0.99851920367213542, 1.0116345126548261, 1.0335801697615423, 1.0116345126548261, 0.99851920367213542, 0.98934385044409012, 0.98371234834817733},
+  optimal_lscale_10[] = {0.98713874330193863, 0.99110664000610571, 0.99544593877259457, 1.0054188667757713, 1.02088981114359, 1.02088981114359, 1.0054188667757713, 0.99544593877259457, 0.99110664000610571, 0.98713874330193863},
+  optimal_lscale_11[] = {0.98802560946845452, 0.99088308151773186, 0.99485717770327386, 1.0018012691703861, 1.0132746372828709, 1.022316449714566, 1.0132746372828709, 1.0018012691703861, 0.99485717770327386, 0.99088308151773186, 0.98802560946845452},
+  optimal_lscale_12[] = {0.98933384410274483, 0.99148492407376065, 0.99456924585078266, 0.99974133924938646, 1.0076339300810422, 1.0172367166422831, 1.0172367166422831, 1.0076339300810422, 0.99974133924938646, 0.99456924585078266, 0.99148492407376065, 0.98933384410274483},
+  optimal_lscale_13[] = {0.99268113885731557, 0.99445091830329768, 0.99613360268945739, 0.99536834490944437, 1.0054811657260543, 1.0081328020609224, 1.0155040549070158, 1.0081328020609224, 1.0054811657260543, 0.99536834490944437, 0.99613360268945739, 0.99445091830329768, 0.99268113885731557},
+  optimal_lscale_14[] = {0.99091904296866629, 0.99240628966641609, 0.99441004182561787, 0.99649973126893876, 1.0006824544025643, 1.0087390373686946, 1.016343402499102, 1.016343402499102, 1.0087390373686946, 1.0006824544025643, 0.99649973126893876, 0.99441004182561787, 0.99240628966641609, 0.99091904296866629};
+
+#endif

src/jacobi/impls/tensor/tensor.c

@@ -11,6 +11,7 @@
 
 #include "tensor.h"
 #include "inlinepoly.h"
+#include "optimalscale.h"
 #include "dohpgeom.h"
 #include "dohpmesh.h"           /* for iMesh_TopologyName */
 
@@ -34,6 +35,21 @@ static dErr TensorJacobiHasBasis(dJacobi,dInt,dInt,dBool*);
 
 static dErr dRealTableView(dInt m,dInt n,const dReal mat[],const char name[],dViewer viewer);
 
+static dErr dJacobiSetFromOptions_Tensor(dJacobi jac)
+{
+  Tensor this = jac->impl;
+  dErr err;
+
+  dFunctionBegin;
+  err = PetscOptionsHead("Tensor options");dCHK(err);
+  {
+    err = PetscOptionsTruth("-djac_tensor_mass_scale","Use optimal scaling for Q1 mass matrices","Jacobi",this->usemscale,&this->usemscale,NULL);dCHK(err);
+    err = PetscOptionsTruth("-djac_tensor_laplace_scale","Use optimal scaling for Q1 Laplacian matrices","Jacobi",this->uselscale,&this->uselscale,NULL);dCHK(err);
+  }
+  err = PetscOptionsTail();dCHK(err);
+  dFunctionReturn(0);
+}
+
 /**
 * Prepare the dJacobi context to return dRule and dEFS objects (with dJacobiGetRule and dJacobiGetEFS).
 *
@@ -71,6 +87,8 @@ static dErr dJacobiSetUp_Tensor(dJacobi jac)
         continue;
       }
       err = TensorBasisCreate(this->basisBuilder,rule,j,&this->basis[i*N+j]);dCHK(err);
+      if (!this->usemscale) this->basis[i*N+j]->mscale = optimal_ones;
+      if (!this->uselscale) this->basis[i*N+j]->lscale = optimal_ones;
     }
   }
   err = dJacobiRuleOpsSetUp_Tensor(jac);dCHK(err);
@@ -666,6 +684,27 @@ static dErr TensorBasisCreate(TensorBuilder build,const TensorRule rule,dInt P,T
       b->derivTranspose[j*Q+i] = b->deriv[i*P+j];
     }
   }
+  switch (P) {
+#define _C(p) case p: b->mscale = optimal_mscale_ ## p; b->lscale = optimal_lscale_ ## p; break
+    _C(2);
+    _C(3);
+    _C(4);
+    _C(5);
+    _C(6);
+    _C(7);
+    _C(8);
+    _C(9);
+    _C(10);
+    _C(11);
+    _C(12);
+    _C(13);
+    _C(14);
+#undef _C
+    default:
+      b->mscale = optimal_ones;
+      b->lscale = optimal_ones;
+      dERROR(1,"optimal scaling not available for this order, this should just be a PetscInfo warning");
+  }
   *basis = b;
   dFunctionReturn(0);
 }
@@ -750,7 +789,7 @@ dErr dJacobiCreate_Tensor(dJacobi jac)
 {
   static const struct _dJacobiOps myops = {
     .SetUp              = dJacobiSetUp_Tensor,
-    .SetFromOptions     = 0,
+    .SetFromOptions     = dJacobiSetFromOptions_Tensor,
     .Destroy            = dJacobiDestroy_Tensor,
     .View               = dJacobiView_Tensor,
     .PropogateDown      = dJacobiPropogateDown_Tensor,
@@ -762,22 +801,28 @@ dErr dJacobiCreate_Tensor(dJacobi jac)
   };
   TensorRuleOptions ropt;
   TensorBasisOptions bopt;
+  Tensor tensor;
   dErr err;
 
   dFunctionBegin;
   err = dMemcpy(jac->ops,&myops,sizeof(struct _dJacobiOps));dCHK(err);
-  err = dNew(struct s_Tensor,&jac->impl);dCHK(err);
+  err = dNew(struct s_Tensor,&tensor);dCHK(err);
   err = dNew(struct s_TensorRuleOptions,&ropt);dCHK(err);
   err = dNew(struct s_TensorBasisOptions,&bopt);dCHK(err);
 
+  tensor->usemscale = dFALSE;
+  tensor->uselscale = dFALSE;
+
   ropt->alpha  = 0.0;
   ropt->beta   = 0.0;
   ropt->family = GAUSS;
-  ((Tensor)jac->impl)->ruleOpts = ropt;
+  tensor->ruleOpts = ropt;
 
   bopt->alpha  = 0.0;
   bopt->beta   = 0.0;
   bopt->family = GAUSS_LOBATTO;
-  ((Tensor)jac->impl)->basisOpts = bopt;
+  tensor->basisOpts = bopt;
+
+  jac->impl = tensor;
   dFunctionReturn(0);
 }

src/jacobi/impls/tensor/tensor.h

@@ -64,6 +64,7 @@ struct s_TensorBasis {
   dInt  P,Q;
   dReal *interp,*deriv,*node;
   dReal *interpTranspose,*derivTranspose;
+  const dReal *mscale,*lscale;        /**< Used to produce optimal scaling of sparse mass and Laplacian matrices */
 };
 
 typedef struct s_TensorBasisOptions *TensorBasisOptions;
@@ -103,6 +104,7 @@ struct s_Tensor {
                                   * then \f$ f(q_i) = \sum_{j=0}^n B[i*n+j] f(x_j) \f$ */
   dInt M;                       /**< number of rules */
   dInt N;                       /**< basis size limit */
+  dTruth usemscale,uselscale;
   /* dBufferList data; */
   struct _dRuleOps *ruleOpsLine,*ruleOpsQuad,*ruleOpsHex;
   struct _dEFSOps *efsOpsLine,*efsOpsQuad,*efsOpsHex;