NAG Technical Report 3/2009

Calling NAG Library Routines from Octave

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).

Contents

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 nag_cmplx_lineqn.cc:

   #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;
         INIT_FAIL(fail);
         // Assign the appropriate value of trans
         if (transpose.compare("Nag_NoTrans") == 0) trans=Nag_NoTrans;
         else if (transpose.compare("Nag_Trans") == 0) trans=Nag_Trans;
         else if (transpose.compare("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
         nag_zgetrf(Nag_RowMajor,n,n,a_nag,pda,ipiv,&fail);
         
         nag_zgetrs(Nag_RowMajor,trans,n,nrhs,a_nag,pda,ipiv,b_nag,pdb,&fail);
         
         // 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:

   • For clarity, the input argument type checking and error message printing have been omitted in the article text. However, they have been added to the file you can download from this site.

   • Similarly to the previous example, we are using the Complex type which is defined in both Octave and NAG header files. Again, we rename it to OctComplex and NagComplex and copy the values from one type to another. This time however we need to copy the whole arrays and matrices of Complex type. We are not using MatrixType, so it can be renamed in one place only, like in the previous examples.

   • nag_stdlib.h needs to be included for the declaration of NAG_ALLOC and NAG_FREE, which select suitable memory allocation and freeing functions.

   • The third DEFUN_DLD argument (nargout) is not used, so it is omitted in order to avoid the warning from gcc about an unused function parameter.

   • #define A(I,J) a_nag[I*tda + J] and #define B(I,J) b_nag[I*tdb + J] are macros to do subscripting into 2 dimensional arrays a_nag and b_nag, which are passed as double a_nag[] and double b_nag[].

3. Compiling into Oct-File

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

     % mkoctfile nag_cmplx_lineqn.cc -L/opt/NAG/cll6a08dg/lib -lnagc_nag -I/opt/NAG/cll6a08dg/include
   

   where:
   
   • the -L switch tells the C compiler where to look for NAG C Library installed on your system;
     
   • -lnagc_nag is the name of the NAG C Library;
     
   • the -I switch should be followed by the location of NAG C Library header files.

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).