NAG Technical Report 3/2009

Calling NAG Library Routines from Octave

Start of report   Skip to Example 1   Skip to Example 2   Skip to Example 4   Skip to Example 5

3.3. Example 3


nag_zgetrs (f07asc) routine - solving a complex system of linear equations AX=B, ATX=B or AHX=B, where A has been factorized by nag_zgetrf (f07arc).

We show how to call a function with complex input and output arrays.


  1. Function prototype from the NAG C Library Manual
  2. C++ function
  3. Compiling into Oct-File
  4. Calling the function
  1. Function prototype from the NAG C Library Manual

    According to the C Library Manual, the prototype for routine f07arc looks like this:

      #include <nag.h>
      #include <nagf07.h>
      void nag_zgetrf(Nag_OrderType order, Integer m, Integer n, Complex a[], Integer pda, Integer ipiv[], NagError *fail);
    And for f07asc:
      #include <nag.h>
      #include <nagf07.h>
      void nag_zgetrs(Nag_OrderType order, Nag_TransType trans, Integer n, Integer nrhs, const Complex a[], Integer pda, const Integer ipiv[], Complex b[], Integer pdb, NagError *fail);
    For argument descriptions consult f07arc and f07asc documents in the NAG C Library Manual [6]. To keep things simple, we are going to call f07arc from within our C++ function. This means that we will not need to pass the ipiv array to and from the function.
    Argument trans is of type Nag_TransType, which does not exist in Octave, so we will pass it as a string argument and then use that string to assign the appropriate value of trans.
    Argument order specifies the two-dimensional storage scheme being used, i.e. row-major ordering or column-major ordering. Octave's defined storage is specified by order=Nag_RowMajor. We will declare this variable inside our C++ function. We will also return only the integer value of fail.code.
  2. C++ function

    Here is the source code of our C++ function
    #undef Complex
    #define Complex OctComplex
    #include <octave/oct.h>
    #undef Complex
    #define Complex NagComplex
    #define MatrixType NagMatrixType
    #include <nag.h>
    #include <nag_stdlib.h>
    #include <nagf07.h>
    DEFUN_DLD (nag_cmplx_lineqn, args, ,
       "Calls nag_zgetrf (f07arc) followed by \n\
    nag_zgetrs (f07asc).\n\
    nag_zgetrs solves a complex system of linear equations with multiple \n\
    right-hand sides, AX=B, A^TX=B or A^HX=B, where A has been factorized \n\
    by nag_zgetrf.\n")
      // Variable to store function output values
      octave_value_list retval;
      // Retrieve input arguments from args
      std::string transpose = args(0).string_value ();
      Integer n = args(1).int_value();
      Integer nrhs = args(2).int_value();
      ComplexMatrix a(args(3).complex_matrix_value());
      Integer pda = args(4).int_value();
      ComplexMatrix b(args(5).complex_matrix_value());
      Integer pdb = args(6).int_value();  
      if (! error_state)
          // Declare local variables
    #define A(I,J) a_nag[I*pda + J]
    #define B(I,J) b_nag[I*pdb + J]
          int i, j;
          NagComplex *a_nag=0, *b_nag=0;
          Integer *ipiv=0;
          Nag_TransType trans;
          NagError fail;
          // Assign the appropriate value of trans
          if ("Nag_NoTrans") == 0) trans=Nag_NoTrans;
          else if ("Nag_Trans") == 0) trans=Nag_Trans;
          else if ("Nag_ConjTrans") == 0) trans=Nag_ConjTrans;
          else error ("transpose not Nag_NoTrans, Nag_Trans or Nag_ConjTrans");
          if ( !( ipiv = NAG_ALLOC(n, Integer))  ||
               !( a_nag = NAG_ALLOC(n*pda, NagComplex)) ||
               !( b_nag = NAG_ALLOC(n*pdb, NagComplex)) )
              error ("Allocation failure\n");     
          // Copy the input OctComplex matrices into NagComplex arrays
          for (i = 0; i < n; i++)
            for (j = 0; j < pda; j++)
                A(i,j).re = a.checkelem(i,j).real();
                A(i,j).im = a.checkelem(i,j).imag();
          for (i = 0; i < n; i++)
            for (j = 0; j < pdb; j++)
                B(i,j).re = b.checkelem(i,j).real();
                B(i,j).im = b.checkelem(i,j).imag();
          // Call NAG routines
          // Copy the output NagComplex array back into OctComplex matrix
          for (i = 0; i < n; i++)
            for (j = 0; j < pdb; j++)
                b.checkelem(i,j).real() = B(i,j).re;
                b.checkelem(i,j).imag() = B(i,j).im; 
          if (a_nag) NAG_FREE(a_nag);
          if (b_nag) NAG_FREE(b_nag);
          if (ipiv) NAG_FREE(ipiv);
          // Assign output arguments to retval
          retval(0) = b;
          retval(1) = fail.code;
      return retval;

    Points to note about this code:

  3. Compiling into Oct-File

    To compile the C++ function into oct-files, we use the mkoctfile script supplied with Octave:

      % mkoctfile -L/opt/NAG/cll6a08dg/lib -lnagc_nag -I/opt/NAG/cll6a08dg/include
  4. Calling the function

    Assuming that all has gone well, we can call the function as if it was part of Octave itself, i.e. either from the Octave command line or from within an Octave program. An example call may be:
    octave:1> a = [-1.3400+2.5500i 0.2800+3.1700i -6.3900-2.2000i 0.7200-0.9200i; -0.1700-1.4100i 3.3100-0.1500i -0.1500+1.3400i 1.2900+1.3800i; -3.2900-2.3900i -1.9100+4.4200i -0.1400-1.3500i 1.7200+1.3500i; 2.4100+0.3900i -0.5600+1.4700i -0.8300-0.6900i -1.9600+0.6700i]
    a =
     Columns 1 through 3:
      -1.34000 + 2.55000i   0.28000 + 3.17000i  -6.39000 - 2.20000i
      -0.17000 - 1.41000i   3.31000 - 0.15000i  -0.15000 + 1.34000i
      -3.29000 - 2.39000i  -1.91000 + 4.42000i  -0.14000 - 1.35000i
       2.41000 + 0.39000i  -0.56000 + 1.47000i  -0.83000 - 0.69000i
     Column 4:
       0.72000 - 0.92000i
       1.29000 + 1.38000i
       1.72000 + 1.35000i
      -1.96000 + 0.67000i
    octave:2> b=[26.2600+51.7800i 31.3200-6.7000i; 6.4300-8.6800i 15.8600-1.4200i; -5.7500+25.3100i -2.1500+30.1900i; 1.1600+2.5700i -2.5600+7.5500i]
    b =
       26.2600 + 51.7800i   31.3200 -  6.7000i
        6.4300 -  8.6800i   15.8600 -  1.4200i
       -5.7500 + 25.3100i   -2.1500 + 30.1900i
        1.1600 +  2.5700i   -2.5600 +  7.5500i
    octave:3> [b,ifail]=nag_cmplx_lineqn("Nag_NoTrans",4,2,a,4,b,2)
    b =
       1.0000 + 1.0000i  -1.0000 - 2.0000i
       2.0000 - 3.0000i   5.0000 + 1.0000i
      -4.0000 - 5.0000i  -3.0000 + 4.0000i
       0.0000 + 6.0000i   2.0000 - 3.0000i
    ifail = 0

    (If you get an error message saying that a library cannot be located, see the tip given in Example 1).

Start of report   Skip to Example 1   Skip to Example 2   Skip to Example 4   Skip to Example 5