2

I am using solvers of LAPACK libraries in a MATLAB MEX file for solving linear system of equations. For some of the cases, the system that I solve is singular. For example, the system is as follows for one of my cases:

A =
 0.00000000 0.00000000  0.00000000
 0.00000000 0.00000000  0.00000000
 0.00000000 0.00000000  77.31867171

b:
-0.00000000 -0.00000000 -0.00000000

What would be the best approach to label the solution of Ax=b of the above system as NaN similar to MATLAB?

iikkoo
  • 2,810
  • 6
  • 33
  • 38
afp_2008
  • 1,940
  • 1
  • 19
  • 46

1 Answers1

4

Here is an example to create a numeric vector filled with NaNs from a MEX-function:

test_nan.cpp

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    plhs[0] = mxCreateDoubleMatrix(3, 1, mxREAL);
    double *x = mxGetPr(plhs[0]);
    double nanVal = mxGetNaN();
    for (int i=0; i<3; ++i) {
        x[i] = nanVal;
    }
}

MATLAB

>> mex -largeArrayDims test_nan.cpp
>> x = test_nan()
x =
   NaN
   NaN
   NaN
Amro
  • 123,847
  • 25
  • 243
  • 454
  • Thanks. As a side question, how would this be done if it was a pure `c++` code? – afp_2008 Jun 24 '14 at 20:56
  • 1
    C++ has a [standard function](http://en.cppreference.com/w/cpp/numeric/math/nan) for that. The thing you have to understand is that in the [IEEE 754 floating-point](https://en.wikipedia.org/wiki/IEEE_754-1985#NaN) representation, there are many values that denote NaN (not unique): http://stackoverflow.com/a/20723890/97160 – Amro Jun 24 '14 at 21:02
  • 1
    @A2009 There's also [`std::numeric_limits::quiet_NaN()`](http://en.cppreference.com/w/cpp/types/numeric_limits/quiet_NaN) in ``. – chappjc Jun 25 '14 at 00:08
  • 2
    @chappjc: note the phrase `Only meaningful if std::numeric_limits::has_quiet_NaN == true`, so there is still a chance (really tiny!) that function wont work for some platforms/implementations out there: http://stackoverflow.com/a/16691798/97160 – Amro Jun 25 '14 at 00:26
  • 1
    @Amro While `has_quiet_NaN` does check that your system fully supports the IEEE 754 standard ([chances that it doesn't are really tiny indeed](http://stackoverflow.com/a/2234499/2778484)), I think its primary purpose is to handle integer types. Since `numeric_limits` is a template class, there's nothing stopping a user from requesting `std::numeric_limits::quiet_NaN()`. What could that even mean, right? Yet VC10 has no problem compiling it. At least they had the sense to specialize that case and return 0. Anyway, the `cmath` functions don't have that issue. – chappjc Jun 25 '14 at 00:42
  • 1
    I just noticed that the [return value for `std::nan`](http://en.cppreference.com/w/cpp/numeric/math/nan#Return_value) is "The quiet NaN value ... _or zero if the implementation does not support quiet NaNs._" So, it sounds like they probably do the same thing on the same system. :) – chappjc Jun 25 '14 at 00:52