topic Does "Sometimes it doesn't in Intel® oneAPI Math Kernel Library

djacobi doesn't exit

Gianluca_G_1 — Wed, 05 Apr 2017 13:41:48 GMT

(Nonlinear Least Squares Problem without Constraints)
I have implemented a C# TrustRegion Algotithm using MKL API (eg. dtrnlspbc_solve, djacobi ...) for our purposes.
It works very well but sometimes no.

Sometimes it doesn't exit, it seems djacobi remains in loop (at least it seems debugging).
In the attached image you can see the CPU usage, I've executed the method three times.
The first time the method terminate, the CPU usage is restored to the previous value,
and so the second time.
But the third time something remains in usage!

I have created a dedicated class with a method that execute the algorithm as a Thread.
After 60 seconds, if the solver hadn't finished it kills the thread, but it doesn't work.

I've tried different tricks but anyone of this solved the problem.

Any suggestion?

Thank you very much

Below the solver main loop ...

while (!bSuccessful)
            {
                if (mbShowMessages) Console.WriteLine("Step " + iStep.ToString());

                if (RCI_Solve(ref handle, fvec, fjac, ref iRCI_Request) != TR_SUCCESS)
                {
                    Console.WriteLine("RCI_Solve() retun error!");
                    RCI_FreeBuffers();
                    iError = 1;
                    goto end;
                }
                bSuccessful =
                    iRCI_Request == -1 || iRCI_Request == -2 || iRCI_Request == -3 ||
                    iRCI_Request == -4 || iRCI_Request == -5 || iRCI_Request == -6;

                if (iRCI_Request == 1)
                {
                    if (mbShowMessages) Console.WriteLine("go objective_function()");
                    objective_function(ref m, ref n, x, fvec);
                }
                if (iRCI_Request == 2)
                {
                    if (mbShowMessages) Console.WriteLine("go djacobi()");
                    int iRes = djacobi(objective_function, ref n, ref m, fjac, x, ref jac_eps);
                    if (iRes != TR_SUCCESS)
                    {
                        if (iRes == TR_INVALID_OPTION)
                        {
                            Console.WriteLine("error in djacobi: invalid options.");
                        }
                        else if (iRes == TR_OUT_OF_MEMORY) 
                        {
                            Console.WriteLine("error in djacobi: out of memory.");
                        }
                       
                        RCI_FreeBuffers();
                        iError = 1;
                        goto end;
                    }

                }

                if (iStep > MAX_STEP)
                {
                    Console.WriteLine(string.Format("Too many external loop! (max {0})", MAX_STEP));
                    RCI_FreeBuffers();
                    iError = 1;
                    goto end;
                }

                iStep++;
            }

Does "Sometimes it doesn't

Jing_Xu — Fri, 07 Apr 2017 08:13:43 GMT

Does "Sometimes it doesn't exit" means that djacobi didn't return? Or it returned but CPU usage was still high?

Yes it didn't return. In this

Gianluca_G_1 — Fri, 07 Apr 2017 08:57:05 GMT

Yes it didn't return. In this case I kill the thread, but it seems it is still working.

Gianluca

It seems that every time I

Gianluca_G_1 — Fri, 07 Apr 2017 13:32:54 GMT

It seems that every time I execute the procedure the cpu time increases.

After three times it doesn't return.

any idea?

How do you call MKL functions

Jing_Xu — Tue, 11 Apr 2017 00:54:41 GMT

How do you call MKL functions in C#?

A workaround is calling these MKL functions inside a cpp file, then build it to be a dll. In C#, you may call the functions implemented in the dll to avoid calling MKL functions directly in C#.

It is long time I don't work

Gianluca_G_1 — Tue, 11 Apr 2017 10:09:45 GMT

It is long time I don't work with c++, that it will take me too many time.

But do you know this problem?

If I share with you a c# console program example, could you help me?

We can use an ufficial support channel, we have a regular licence of Parallel studio XE.

Gianluca

Sometimes, it is better to

mecej4 — Tue, 11 Apr 2017 11:59:01 GMT

Sometimes, it is better to break up the problem into smaller pieces. Let us set aside the complexities of mixed language calls between C#, C++ and MKL. What is the nature of this objective function for which it takes ~60 seconds to compute the Jacobian? How many functions make up the SSQ expression, how many variables, and what is the time for one SSQ evaluation?

Image a button on a Form that

Gianluca_G_1 — Tue, 11 Apr 2017 12:19:56 GMT

Image a button on a Form that I click every time I need the calculation.

The question is that the method works well the first time I call it. Also the second time but the third not.

It doesn't depend on the complexity of "obj func" at least it seems like this.

Below L3_Function(), one case of objective function with 5 variables.

It gives a Soil Model Resistivity aproximation. If you please I can be more precise.

public void L3_Function(ref int m, ref int n, [In] IntPtr xp, [Out] IntPtr fp)
        {
            double[] x = new double;
            Marshal.Copy(xp, x, 0, n);
            double[] f = new double;
            Marshal.Copy(fp, f, 0, m);

            double rho1 = x[0];
            double rho2 = x[1];
            double rho3 = x[2];
            double h1 = x[3];
            double h2 = x[4];
            double w = 1;

            for (int i = 0; i < m; i++)
            {
                w = Math.Sqrt(mlMeasures.NumA);
                if (mbApplyMWeight)
                {
                    w = Math.Sqrt((1.0 / m + Math.Abs(i / m - 1 / 2)) * mlMeasures.NumA);
                }

                f = (1 - L3_RhoCalc(mlMeasures.A, rho1, rho2, rho3, h1, h2) / mlMeasures.RhoM) * w;
            }

            Marshal.Copy(f, 0, fp, m);
        }

        public double L3_RhoCalc(double a, double rho1, double rho2, double rho3, double h1, double h2)
        {          
            double dRe2 = AdpGauss.ThreeLayers(a, rho1, rho2, rho3, h1, h2);
            return dRe2;
        }

public static double ThreeLayers(double a, double rho1, double rho2, double rho3, double h1, double h2)
        {
            double dError = MAX_ERROR;
            int iStep = STEP_BASE;
            int iCount = 0;
            double dLeft = 0;
            double dRight = 0;
            double dInt = 0; //Integral Value
            double dIntPrev = 0; //Integral Previous Value

            double dNi12 = (rho2 - rho1) / (rho2 + rho1);
            double dNi23 = (rho3 - rho2) / (rho3 + rho2);

            double[] w = new double[] { 5d / 9d, 8d / 9d, 5d / 9d };
            double[] x = new double[] { -1d * Math.Sqrt(3d / 5d), 0d, Math.Sqrt(3d / 5d) };

            #region CalcLeft
            while (dError > EPSILON)
            {
                if (iCount > MAX_ITERATIONS)
                    break;

                dLeft = 0;
                //Console.WriteLine("> Step = " + iStep.ToString());
                for (int i = 0; i < iStep; i++)
                {
                    for (int j = 0; j < 3; j++)
                    {
                        double t = (x + 2 * i + 1) / (2 * iStep);
                        //Console.WriteLine(string.Format("t({0},{1}) =  {2}", i, j, t.ToString()));
                        double k31 = (dNi12 + dNi23 * Math.Pow(Math.E, -2 * t * h2)) / (1 + dNi12 * dNi23 * Math.Pow(Math.E, -2 * t * h2));
                        double b3 = k31 * Math.Pow(Math.E, -2 * t * h1) / (1 - k31 * Math.Pow(Math.E, -2 * t * h1));
                        double f = b3 * (RciWrapper.Bessel(t * a) - RciWrapper.Bessel(2 * t * a));
                        dLeft += f * w;
                    }
                }

                dLeft /= 2 * iStep;
                dInt = dLeft;

                if (iCount > 0)
                {
                    dError = Math.Abs(1 - dIntPrev / dInt);
                    //Console.WriteLine(iCount.ToString() + "> AdpGauss Error: " + dError.ToString());
                }

                dIntPrev = dInt;
                iStep *= 2;
                iCount++;
            }
            #endregion

            #region CalcRight
            dError = MAX_ERROR;
            iStep = 1;
            iCount = 0;
            while (dError > EPSILON)
            {
                if (iCount > MAX_ITERATIONS)
                    break;

                dRight = 0;
                for (int i = 0; i < iStep; i++)
                {
                    for (int j = 0; j < 3; j++)
                    {
                        double t = (x + 2 * i + 1) / (2 * iStep);
                        double k31 = (dNi12 + dNi23 * Math.Pow(Math.E, -2 * (1 / t) * h2)) / (1 + dNi12 * dNi23 * Math.Pow(Math.E, -2 * (1 / t) * h2));
                        double b3 = k31 * Math.Pow(Math.E, -2 * (1 / t) * h1) / (1 - k31 * Math.Pow(Math.E, -2 * (1 / t) * h1));
                        double f = b3 * (RciWrapper.Bessel((1 / t) * a) - RciWrapper.Bessel(2 * (1 / t) * a));
                        dRight += (f * w) / (t * t);
                    }
                }

                dRight /= 2 * iStep;

                dInt = dRight;

                if (iCount > 0)
                {
                    dError = Math.Abs(1 - dIntPrev / dInt);
                    //Console.WriteLine(iCount.ToString() + "> AdpGauss Error: " + dError.ToString());
                }

                dIntPrev = dInt;
                iStep *= 2;
                iCount++;
            }
            #endregion

            return rho1 + 4 * rho1 * a * (dLeft + dRight);
        }

Could you provide your

Jing_Xu — Wed, 12 Apr 2017 02:11:40 GMT

Could you provide your complete code, please? Or a complete sample code.

Is it possible to send it in

Gianluca_G_1 — Wed, 12 Apr 2017 10:21:45 GMT

Is it possible to send it in a confidential way?