Method DiffJacobi

fem/coefficient_matrix.cpp:435–459

    }

    shared_ptr<CoefficientFunction> DiffJacobi(const CoefficientFunction *var,
                                               T_DJC &cache) const override
    {
      auto thisptr = const_pointer_cast<CoefficientFunction>(this->shared_from_this());
      if (cache.find(thisptr) != cache.end())
        return cache[thisptr];

      if (this == var)
        return IdentityCF(this->Dimensions());

      auto diffc1 = c1->DiffJacobi (var, cache);
      auto inv1 = thisptr;

      const int D = c1->Dimensions()[0];
      Array<int> dimres { D, D };
      dimres += var->Dimensions();

      auto prod1 = -inv1 * diffc1->Reshape( D, -1 );
      auto prod1r = prod1 -> Reshape( Array<int> (dimres) );
      auto trans = prod1r -> TensorTranspose( 0, 1 );
      auto prod2 = inv1->Transpose() * trans -> Reshape( D, -1 );
      auto res = prod2 -> Reshape( Array<int> (dimres) ) -> TensorTranspose( 0, 1);
      cache[thisptr] = res;
      return res;
    }
  };