C++ DoubleMatrix::clear方法代码示例

double IFunctional::fx(const DoubleVector &pv) const
{
    IFunctional* ifunc = const_cast<IFunctional*>(this);

    ifunc->fromVector(pv, ifunc->mParameter);
    //functional->mParameter.fromVector(prms);
    ifunc->forward->setParameter(mParameter);
    ifunc->backward->setParameter(mParameter);

    DoubleMatrix u;
    vector<ExtendedSpaceNode2D> info;
    ifunc->forward->calculateMVD(u, info, true);

    double intgrl = integral(u);

    u.clear();

    return intgrl + regEpsilon*norm() + r*penalty(info);
}

void BackwardParabolicIBVP::gridMethod(DoubleMatrix &p, SweepMethodDirection direction) const
{
    typedef void (*t_algorithm)(const double*, const double*, const double*, const double*, double*, unsigned int);
    t_algorithm algorithm = &tomasAlgorithm;
    if (direction == ForwardSweep) algorithm = &tomasAlgorithmL2R;
    if (direction == BackwardSweep) algorithm = &tomasAlgorithmR2L;

    Dimension time = mtimeDimension;
    Dimension dim1 = mspaceDimension.at(0);

    double ht = time.step();
    unsigned int minM = time.min();
    unsigned int maxM = time.max();
    unsigned int M = maxM-minM;

    double hx = dim1.step();
    unsigned int minN = dim1.min();
    unsigned int maxN = dim1.max();
    unsigned int N = maxN-minN;

    double h = ht/(hx*hx);

    p.clear();
    p.resize(M+1, N+1);

    double *ka = (double*) malloc(sizeof(double)*(N-1));
    double *kb = (double*) malloc(sizeof(double)*(N-1));
    double *kc = (double*) malloc(sizeof(double)*(N-1));
    double *kd = (double*) malloc(sizeof(double)*(N-1));
    double *rx = (double*) malloc(sizeof(double)*(N-1));

    /* initial condition */
    SpaceNodePDE isn;
    for (unsigned int n=0; n<=N; n++)
    {
        isn.i = n+minN;
        isn.x = isn.i*hx;
        p[M][n] = initial(isn);
    }
    layerInfo(p, M);

    SpaceNodePDE lsn;
    lsn.i = minN;
    lsn.x = minN*hx;

    SpaceNodePDE rsn;
    rsn.i = maxN;
    rsn.x = maxN*hx;

    TimeNodePDE tn;
    for (unsigned int m=M-1; m!=UINT32_MAX; m--)
    {
        tn.i = m+minM;
        tn.t = tn.i*ht;

        for (unsigned int n=1; n<=N-1; n++)
        {
            isn.i = n+minN;
            isn.x = isn.i*hx;

            double alpha = -a(isn,tn)*h;
            double betta = 1.0 - 2.0*alpha;

            ka[n-1] = alpha;
            kb[n-1] = betta;
            kc[n-1] = alpha;
            kd[n-1] = p[m+1][n] - ht * f(isn, tn);
        }

        ka[0]   = 0.0;
        kc[N-2] = 0.0;

        /* border conditions */
        p[m][0] = boundary(lsn, tn);
        p[m][N] = boundary(rsn, tn);

        kd[0]   += a(lsn,tn) * h * p[m][0];
        kd[N-2] += a(rsn,tn) * h * p[m][N];

        (*algorithm)(ka, kb, kc, kd, rx, N-1);

        for (unsigned int n=1; n<=N-1; n++) p[m][n] = rx[n-1];

        layerInfo(p, m);
    }

    free(ka);
    free(kb);
    free(kc);
    free(kd);
    free(rx);
}

void DiscreteHeat::calculateP(const DoubleVector &f, const DoubleVector &u, DoubleMatrix &psi, DoubleVector &g)
{
    C_UNUSED(f);
    C_UNUSED(g);

    double lamda = -(a*ht)/(hx*hx);
    double k = 1.0-2.0*lamda;

    psi.clear();
    psi.resize(M+1, N+1);

    DoubleVector a1(N-1);
    DoubleVector b1(N-1);
    DoubleVector c1(N-1);
    DoubleVector d1(N-1);
    DoubleVector x1(N-1);

    for (unsigned int m=0; m<=M; m++)
    {
        unsigned int j = M-m;

        if (j==M)
        {
            for (unsigned int i=1; i<=N-1; i++)
            {
                a1[i-1] = lamda;
                b1[i-1] = k;
                c1[i-1] = lamda;
                d1[i-1] = -2.0*hx*(u[i]-U[i]);
            }
            a1[0]   = 0.0;
            c1[N-2] = 0.0;
            tomasAlgorithm(a1.data(), b1.data(), c1.data(), d1.data(), x1.data(), x1.size());

            for (unsigned int i=1; i<=N-1; i++)
            {
                psi[j][i] = x1[i-1];
            }
            psi[j][0] = -lamda*psi[j][1]   -2.0*hx*0.5*(u[0]-U[0]);
            psi[j][N] = -lamda*psi[j][N-1] -2.0*hx*0.5*(u[N]-U[N]);
        }
        else if (j==0)
        {
            psi[j][0] = -lamda*psi[j][1];
            psi[j][N] = -lamda*psi[j][N-1];

            for (unsigned int i=1; i<=N-1; i++)
            {
                psi[j][i] = psi[j+1][i];
            }
        }
        else
        {
            for (unsigned int i=1; i<=N-1; i++)
            {
                a1[i-1] = lamda;
                b1[i-1] = k;
                c1[i-1] = lamda;
                d1[i-1] = +psi[j+1][i];
            }
            a1[0]   = 0.0;
            c1[N-2] = 0.0;
            tomasAlgorithm(a1.data(), b1.data(), c1.data(), d1.data(), x1.data(), x1.size());

            for (unsigned int i=1; i<=N-1; i++)
            {
                psi[j][i] = x1[i-1];
            }
            psi[j][0] = -lamda*psi[j][1];
            psi[j][N] = -lamda*psi[j][N-1];
        }
    }
    a1.clear();
    b1.clear();
    c1.clear();
    d1.clear();
    x1.clear();
}