Using the C++ Math Library

The C++ roots() Function

The example is divided into two parts. The first part shows the main program, which sets up the problem and invokes the example version of roots(): the example_roots() function. The second part contains the example_roots() function. In the C++ source file, the order of the parts is reversed. The parts are reordered here for clarity.

You can find the code for this example in the
<matlab>/extern/examples/cppmath directory on UNIX systems and in the <matlab>\extern\examples\cppmath directory on PCs, where <matlab> represents the top-level directory of your installation. See Building C++ Applications for information about building and running the example program.

In the example, note the following:

• Programs written using the MATLAB C++ Math Library look very much like MATLAB M-files. The syntax of the two is very similar.
• The MATLAB C++ Math Library supports most of MATLAB's operators. Those that are not supported as operators can be accessed via function calls.
• The functions colon() and ramp() in the MATLAB C++ Math Library replace the MATLAB colon operator.
• A function that modifies its input arguments, or uses them as a temporary variables, must not declare those arguments const.
• The default mwIndex constructor produces an mwIndex object that acts, when used as a subscript, like a colon in MATLAB.
• The C++ if-statement, unlike the MATLAB if statement, requires that any matrices tested be reduced in rank to scalars by calls to any() or all(). See the online MATLAB C++ Math Library Reference for further explanation of the functions. Accessing Online Reference Documentation describes how to access the Help Desk.
• Calling the error() function throws an mwRuntimeException exception. mwRuntimeException is a subclass of mwException.
• vertcat() vertically concatenates two or more matrices. horzcat() behaves like vertcat() except performs horizontal concatenation.
• Use the C++ operator ! to invert the truth value of a logical scalar mwArray; use the MATLAB C++ Math Library operator ~ to invert the truth value of a logical array. Do not apply ! to arrays.

// Call example_roots() and the library's roots().
int main(void)
{
    // Static array of doubles used to initialize the matrices.
    static double input[] = { 1, -6, -72, -27 };            // #1

    // Declare three matrices, one with initial values.
    mwArray x(1, 4, input), result, verify;                 // #2

    // Call our version of roots().
    result = example_roots(x);                              // #3

    // Call the MATLAB C++ Math Library roots.
    verify = roots(x);                                      // #4

    // Print the input and output matrices from example_roots().
    cout << "x = " << endl << x << endl;
    cout << "example_roots(x) = " << endl << result << endl; // #5

    // Check to see if the answer is equal to the real roots().
    if (tobool(result == verify))                           // #6
        cout << "Success!" << endl;

    return(EXIT_SUCCESS);
}

The numbered items in the list below correspond to the numbered comments in the code example:

The main program is straightforward, so the explanations below are brief. If you have difficulty understanding this section of the example, refer to Example Program: Writing Simple Functions (ex4.cpp) in Chapter 2.

1. Declare the static variable used for array initialization. The elements of the C++ array are specified in the column-major order required by the library.
2. Declare and initialize three matrices. x, a row vector (one row, four columns), is the input matrix. result and verify are initially null matrices.
3. Call the example_roots() function and place the return value in the matrix result.
4. Call the MATLAB C++ Math Library's version of roots() and store the return value in the matrix verify.  verify will be used to confirm that the rewriting of roots produces the correct result.
5. Print the input to example_roots(). Print the output from example_roots().
6. Verify the result. The matrix verify contains the correct result. The mwArray class provides an overloaded operator==(), which makes comparing two matrices for equality easy.

The second part of the example is the example_roots() function itself. This function is part of the same file as the main program shown above.

#include <stdlib.h>
#include "matlab.hpp"                                      // #1

// EXAMPLE_ROOTS(C) computes the roots of the polynomial whose
// coefficients are the elements of the vector C. If C has N+1
// components, the polynomial is C(1)*X^N + ... + C(N)*X + C(N+1).
mwArray example_roots(mwArray c)                          // #2
{
    mwArray n, inz, nnz, r, a;                            // #3
    mwIndex icolon;

    // Make sure number of dimensions is not greater than 1.
    n = size(c);                                          // #4
    if (all(n > 1.0))
        error("Must be a vector");

    n = max(n);                                          // #5
    c = transpose(c(icolon)); // Make sure it's a row vector. 
    inz = find(abs(c));       // Find nonzero elements.     #6
    nnz = max(size(inz));     // Count nonzero elements.

    // Test all elements against zero.
    if (!(nnz == 0.0))                                   // #7
    {
        c = c(ramp(inz(1), inz(nnz))); //Strip leading/trailing 0's
        r = zeros(n - inz(nnz), 1);    //Remember trailing 0's
    }   

    //  Polynomial roots via a companion matrix
    n = max(size(c));  // Size of largest dimension of c   #8
    a = diag(ones(1, n - 2.0), -1.0); // Create a row vector of 1's.
    if (n > 1.0)                                        // #9
        a(1,icolon) = -c(ramp(2, n)) / c(1);
    r = vertcat(r, eig(a));                             // #10

    return r;
}

The numbered items in the list below correspond to the numbered comments in the code example:

1. Include header files. matlab.hpp declares the MATLAB C++ Math Library's data types and functions. matlab.hpp includes iostream.h, which declares the input and output streams cin and cout.  stdlib.h contains the definition of EXIT_SUCCESS.
2. Declare the C++ function example_roots(). We can't use the name roots() because the MATLAB C++ Math Library defines a roots() function with exactly the same number and type of input and output arguments. example_roots() has one input and one output, both of which are matrices. The input argument is not declared const because example_roots() stores temporary results in c.
3. Declare the matrix and index variables. All variables used in the program must be declared before being used. Since icolon is declared using the default mwIndex constructor, this variable acts like MATLAB's : operator when used in array indexing expressions. Declaring a variable like this is more efficient than repeatedly calling the colon() function. See Programming Efficient Indices in Chapter 4 for more information.
4. Check for valid inputs. Determine the size of the input matrix and store the result (a 1-by-2 matrix, i.e., a vector) in the variable n. At least one of the dimensions of the input matrix must be equal to one, that is, the matrix must be a vector. Report an error to the user and terminate if the matrix is malformed. Calling the error() function causes a runtime exception to be thrown.
5. Determine the size of the largest dimension of the input matrix. n is now a 1-by-1 matrix: a scalar. Use the colon operator (here represented by the variable icolon) to extract the input matrix into a column vector. Transpose this vector to get a row vector. The root-finding algorithm below requires that the matrix be a row vector. You could improve the efficiency here by testing the dimensions of the input matrix and transforming them only when necessary.
6. Find all the nonzero elements of the input matrix. Store the result in a vector. Then count the number of nonzero elements. Since inz is a vector, size() returns a vector [ 1 N ], where N is the length of the vector. N is the count of elements in inz.
7. Strip leading zeros and delete them. Strip trailing zeros, but remember them as roots at zero. It is possible that the input matrix was full of zeros. In this case, find() will have returned a null matrix and nnz will be equal to 0. Note that wrapping the logical expression with all() is not necessary in this case, since nnz is known to be a scalar.

If the input matrix did not contain all zeros, inz contains the nonzero elements. Replace the input matrix with a vector 1,2,...,N where N is the number of nonzero elements originally in the input matrix. The result from the call to ramp() goes from the first nonzero index to the last.

Set r to a column vector of zeros, with one row for each trailing zero element of the input matrix. The arguments to the zeros() function are row count and column count.

8. Determine the size of the largest dimension of the input matrix, which may have changed, since the zero elements have been removed from the matrix. Use this size to form a diagonal matrix, with 1's on the -1 diagonal.
9. If c is not empty, replace the first row of matrix a with the vector resulting from dividing the negative of the 2,...,N elements of c (enumerated by the call to ramp()) by c(1). This type of assignment, where an indexing expression appears on the left-hand side, is the only way to modify the contents of a matrix.
10. Vertically concatenate the matrices r and eig(a). The number of columns in both matrices must be the same. The rows of eig(a) are placed below the rows (in this case, the single row) of r. Reassign the result to r. Return the matrix r. Unlike MATLAB, C++ requires an explicit return statement.

Output

The program produces this output:

x = 
  [
         1       -6      -72      -27
  ]
 
 
example_roots(x) = 
  1.0e+01 *
 
  [
    1.21229 ;
   -0.57345 ;
   -0.03884
  ]
 
 
Success!

The M-File roots() Function

mwArray Class Interface