Use 0-indexed vectorized AM in LUCIN

   * IER     IS REPLACED BY 1 FOR SINGULARITY.

   */

  double *v = new double[nddmst];

  dcomplex* vec_am = am[0];

  dcomplex ctemp, cfun;

  const dcomplex cc0 = 0.0 + 0.0 * I;

  ier = 0;

  int nminus = n - 1;

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

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

    double sum = 0.0;

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

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

      sum += (

	      real(am[i - 1][j - 1]) * real(am[i - 1][j - 1])

	      + imag(am[i - 1][j - 1]) * imag(am[i - 1][j - 1])

	      real(vec_am[n * i + j]) * real(vec_am[n * i + j])

	      + imag(vec_am[n * i + j]) * imag(vec_am[n * i + j])

	      );

    } // j1319 loop

    v[i - 1] = 1.0 / sum;

    v[i] = 1.0 / sum;

  } // i1309 loop

  // 2. REPLACE AM BY TRIANGULAR MATRICES (L,U) WHERE AM=L*U.

  //    REPLACE L(I,I) BY 1/L(I,I), READY FOR SECTION 4.

  //    (ROW INTERCHANGES TAKE PLACE, AND THE INDICES OF THE PIVOTAL ROWS

  //    ARE PLACED IN V.)

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

    int64_t kplus = k + 1;

    int64_t kminus = k - 1;

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

    int64_t kplus = k + 2;

    int64_t kminus = k;

    int64_t l = k;

    double psqmax = 0.0;

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

      cfun = cdtp(-am[i - 1][k - 1], am, i, 1, k, kminus);

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

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

      ctemp = -cfun;

      am[i - 1][k - 1] = ctemp;

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

      vec_am[n * i + k] = ctemp;

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

      if (psq > psqmax) {

	psqmax = psq;

	l = i;

      }

    } // i2029 loop

    if (l != k) {

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

	ctemp = am[k - 1][j - 1];

	am[k - 1][j - 1] = am[l - 1][j - 1];

	am[l - 1][j - 1] = ctemp;

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

	ctemp = vec_am[n * k + j];

	vec_am[n * k + j] = vec_am[n * l + j];

	vec_am[n * l + j] = ctemp;

      } // j2049 loop

      v[l - 1] = v[k - 1];

      v[l] = v[k];

    }

    // label 2011

    v[k - 1] = 1.0 * l;

    v[k] = 1.0 * l;

    if (psqmax == 0.0) {

      ier = 1;

      delete[] v;

      return;

    }

    ctemp = 1.0 / am[k - 1][k - 1];

    am[k - 1][k - 1] = ctemp;

    ctemp = 1.0 / vec_am[n * k + k];

    vec_am[n * k + k] = ctemp;

    if (kplus <= n) {

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

	cfun = cdtp(-am[k - 1][j - 1], am, k, 1, j, kminus);

	am[k - 1][j - 1] = -ctemp * cfun;

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

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

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

      } // j2059 loop

    }

  } // k2019 loop

  // 4.  REPLACE AM BY ITS INVERSE AMINV.

  // 4.1 REPLACE L AND U BY THEIR INVERSE LINV AND UINV.

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

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

    int64_t kplus = k + 1;

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

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

      am[i - 1][k - 1] = -am[i - 1][i - 1] * cfun;

      cfun = cdtp(am[k - 1][i - 1], am, k, kplus, i, i - k - 1);

      am[k - 1][i - 1] = -cfun;

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

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

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

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

    } // i4119 loop

  } // k4109 loop

  // 4.2 FORM AMINV=UINV*LINV.

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

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

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

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

      if (i < k) {

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

	am[i - 1][k -1] = cfun;

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

	vec_am[n * i + k] = cfun;

      }

      else {

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

	am[i - 1][k - 1] = cfun;

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

	vec_am[n * i + k] = cfun;

      }

    } // i4119 loop

  } // k4209 loop

  // 4.3 INTERCHANGE COLUMNS OF AMINV AS SPECIFIED BY V, BUT IN REVERSE

  //     ORDER.

  for (int64_t l = 1; l <= n; l++) {

    int64_t k = n - l + 1;

    int64_t kcol = (int64_t)(v[k - 1]);

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

    int64_t k = n - 1 - l;

    int64_t kcol = (int64_t)(v[k]);

    if (kcol != k) {

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

	ctemp = am[i - 1][k - 1];

	am[i - 1][k - 1] = am[i - 1][kcol - 1];

	am[i - 1][kcol - 1] = ctemp;

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

	ctemp = vec_am[n * i + k];

	vec_am[n * i + k] = vec_am[n * i + kcol];

	vec_am[n * i + kcol] = ctemp;

      } // i4319 loop

    }

  } // l4309 loop