Commit f0d4acb7 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Have cdtp() working on the vectorized version of AM

parent 6cea13bc

src/include/clu_subs.h

+2 −1
Original line number Original line Diff line number Diff line
@@ -79,8 +79,9 @@ void apcra( 
 * \param jf: `int`

 * \param jf: `int`

 * \param k: `int`

 * \param k: `int`

 * \param nj: `int`

 * \param nj: `int`

 * \param nddmst: `np_int` Maximum size of the matrix.

 */

 */

dcomplex cdtp(dcomplex z, dcomplex **am, int i, int jf, int k, int nj);

dcomplex cdtp(dcomplex z, dcomplex *vec_am, int i, int jf, int k, int nj, np_int nddmst);





/*! \brief C++ porting of CGEV. QUESTION: description?

/*! \brief C++ porting of CGEV. QUESTION: description?

 *

 *

src/libnptm/clu_subs.cpp

+8 −8
Original line number Original line Diff line number Diff line
@@ -388,7 +388,7 @@ void apcra( 
  delete[] svs;

  delete[] svs;

}

}





dcomplex cdtp(dcomplex z, dcomplex **am, int i, int jf, int k, int nj) {

dcomplex cdtp(dcomplex z, dcomplex *vec_am, int i, int jf, int k, int nj, np_int nddmst) {

  /* NOTE: the original FORTRAN code treats the AM matrix as a

  /* NOTE: the original FORTRAN code treats the AM matrix as a

   * vector. This is not directly allowed in C++ and it requires

   * vector. This is not directly allowed in C++ and it requires

   * accounting for the different dimensions.

   * accounting for the different dimensions.
@@ -397,7 +397,7 @@ dcomplex cdtp(dcomplex z, dcomplex **am, int i, int jf, int k, int nj) { 
  if (nj > 0) {

  if (nj > 0) {

    int jl = jf + nj - 1;

    int jl = jf + nj - 1;

    for (int j = jf; j <= jl; j++) {

    for (int j = jf; j <= jl; j++) {

      result += (am[i - 1][j - 1] * am[j - 1][k - 1]);

      result += (vec_am[nddmst * (i - 1) + j - 1] * vec_am[nddmst * (j - 1) + k - 1]);

    }

    }

  }

  }

  return result;

  return result;
@@ -1181,7 +1181,7 @@ void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier) { 
    int64_t l = k;

    int64_t l = k;

    double psqmax = 0.0;

    double psqmax = 0.0;

    for (int64_t i = k; i < n; i++) {

    for (int64_t i = k; i < n; i++) {

      cfun = cdtp(-vec_am[n * i + k], am, i + 1, 1, k + 1, kminus);

      cfun = cdtp(-vec_am[n * i + k], vec_am, i + 1, 1, k + 1, kminus, n);

      ctemp = -cfun;

      ctemp = -cfun;

      vec_am[n * i + k] = ctemp;

      vec_am[n * i + k] = ctemp;

      double psq = v[i] * (real(ctemp) * real(ctemp) + imag(ctemp) * imag(ctemp));

      double psq = v[i] * (real(ctemp) * real(ctemp) + imag(ctemp) * imag(ctemp));
@@ -1209,7 +1209,7 @@ void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier) { 
    vec_am[n * k + k] = ctemp;

    vec_am[n * k + k] = ctemp;

    if (kplus <= n) {

    if (kplus <= n) {

      for (int64_t j = k + 1; j < n; j++) {

      for (int64_t j = k + 1; j < n; j++) {

	cfun = cdtp(-vec_am[n * k + j], am, k + 1, 1, j + 1, kminus);

	cfun = cdtp(-vec_am[n * k + j], vec_am, k + 1, 1, j + 1, kminus, n);

	vec_am[n * k + j] = -ctemp * cfun;

	vec_am[n * k + j] = -ctemp * cfun;

      } // j2059 loop

      } // j2059 loop

    }

    }
@@ -1219,9 +1219,9 @@ void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier) { 
  for (int64_t k = 0; k < nminus; k++) {

  for (int64_t k = 0; k < nminus; k++) {

    int64_t kplus = k + 1;

    int64_t kplus = k + 1;

    for (int64_t i = kplus; i < n; i++) {

    for (int64_t i = kplus; i < n; i++) {

      cfun = cdtp(cc0, am, i + 1, k + 1, k + 1, i - k);

      cfun = cdtp(cc0, vec_am, i + 1, k + 1, k + 1, i - k, n);

      vec_am[n * i + k] = -vec_am[n * i + i] * cfun;

      vec_am[n * i + k] = -vec_am[n * i + i] * cfun;

      cfun = cdtp(vec_am[n * k + i], am, k + 1, kplus + 1, i + 1, i - k - 1);

      cfun = cdtp(vec_am[n * k + i], vec_am, k + 1, kplus + 1, i + 1, i - k - 1, n);

      vec_am[n * k + i] = -cfun;

      vec_am[n * k + i] = -cfun;

    } // i4119 loop

    } // i4119 loop

  } // k4109 loop

  } // k4109 loop
@@ -1229,11 +1229,11 @@ void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier) { 
  for (int64_t k = 0; k < n; k++) {

  for (int64_t k = 0; k < n; k++) {

    for (int64_t i = 0; i < n; i++) {

    for (int64_t i = 0; i < n; i++) {

      if (i < k) {

      if (i < k) {

	cfun = cdtp(cc0, am, i + 1, k + 1, k + 1, n - k);

	cfun = cdtp(cc0, vec_am, i + 1, k + 1, k + 1, n - k, n);

	vec_am[n * i + k] = cfun;

	vec_am[n * i + k] = cfun;

      }

      }

      else {

      else {

	cfun = cdtp(vec_am[n * i + k], am, i + 1, i + 2, k + 1, n - i - 1);

	cfun = cdtp(vec_am[n * i + k], vec_am, i + 1, i + 2, k + 1, n - i - 1, n);

	vec_am[n * i + k] = cfun;

	vec_am[n * i + k] = cfun;

      }

      }

    } // i4119 loop

    } // i4119 loop