<html>
  <head>

    <meta http-equiv="content-type" content="text/html; charset=UTF-8">
  </head>
  <body>
    <p>I was asked to send out an example of using C++ and Armadillo
      using complex numbers.  I am including here my ztestarma.cpp
      source, which solves the same problem that ztest.f90 and
      gslztest.c solved.  Basically a complex matrix is created, I take
      the inverse of the matrix, and them multiply the inverse and the
      original matrix.  The solution should be a matrix with all 1 along
      the diagonal -- but due to numerical precision errors this is
      rarely the case.  At the end I sum up the diagonal elements to see
      if it matches the matrix dimension.</p>
    <p><font face="monospace">#include <iostream><br>
        #include <armadillo><br>
        // Compile with gcc gslztest.c -lgsl -lm <br>
        <br>
        using namespace std;<br>
        <br>
        int main (void)<br>
        {<br>
        <br>
            const int N = 10;<br>
            int i, j;<br>
            double pi, NCUBE, alpha;<br>
            complex<double> sumdiag;<br>
        <br>
        <br>
            // Create pointers and allocate space for complex matrices<br>
            arma::cx_dmat A(N,N);<br>
            arma::cx_dmat B(N,N);<br>
            arma::cx_dmat C(N,N);<br>
            arma::cx_dmat INV(N,N);<br>
            <br>
        <br>
            pi = acos(-1.0);<br>
        <br>
            // Set the value of N^3 since we use is a lot<br>
        <br>
            NCUBE = (double) N * (double) N * (double) N;<br>
        <br>
            // Now fill the matrices with the same info<br>
        <br>
            for (i=1; i<=N; i++) <br>
                for (j=1; j<=N; j++) {<br>
                    alpha =  sqrt( pow(cos ( (double) (  i*i*j) * pi /
        NCUBE ),2) +<br>
                                   pow(sin ( (double) (2*i*i*j) * pi /
        NCUBE ),2));<br>
                    alpha = (double) 1.0 / alpha; <br>
                    A( i-1, j-1) = complex<double>( <br>
                              alpha * cos( (double) (  i*i*j) * pi /
        NCUBE ),<br>
                              alpha * sin( (double) (2*i*i*j) * pi /
        NCUBE ));<br>
                    B( i-1, j-1) = A( i-1, j-1);<br>
                } <br>
        <br>
            inv(INV,B);<br>
        <br>
            C = INV*A; <br>
           <br>
            sumdiag = complex<double>(0.0,0.0); <br>
            for (i=0;i<N;i++) sumdiag = sumdiag + C(i,i); <br>
        <br>
        <br>
            printf("Sum of diagonal is %g + i %g\n", sumdiag);<br>
        <br>
        }</font><br>
      <br>
    </p>
    <div class="moz-signature">-- <br>
      <b>Andrew J. Pounds, Ph.D.</b><br>
      <i>Professor of Chemistry and Computer Science</i><br>
      <i>Director of the Computational Science Program</i><br>
      <i>Mercer University, Macon, GA, 31207 (478) 301-5627 </i></div>
  </body>
</html>