NAG Library Manual, Mark 30.1
Interfaces:  FL   CL   CPP   AD 

NAG CL Interface Introduction
Example description
/* nag_lapackeig_zggrqf (f08ztc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.1, 2024.
 */

#include <nag.h>
#include <stdio.h>

int main(void) {
  /* Scalars */
  Complex alpha, beta;
  double rnorm;
  Integer i, j, m, n, p, pda, pdb, pdc, pdd, pdx;
  Integer y1rows, y2rows, y3rows;
  Integer exit_status = 0;

  /* Arrays */
  Complex *a = 0, *b = 0, *c = 0, *d = 0, *taua = 0, *taub = 0, *x = 0;

  /* Nag Types */
  NagError fail;
  Nag_OrderType order;

#ifdef NAG_COLUMN_MAJOR
#define A(I, J) a[(J - 1) * pda + I - 1]
#define B(I, J) b[(J - 1) * pdb + I - 1]
  order = Nag_ColMajor;
#else
#define A(I, J) a[(I - 1) * pda + J - 1]
#define B(I, J) b[(I - 1) * pdb + J - 1]
  order = Nag_RowMajor;
#endif

  INIT_FAIL(fail);

  printf("nag_lapackeig_zggrqf (f08ztc) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n]");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &m, &n, &p);
  if (n < 0 || m < 0 || p < 0) {
    printf("Invalid n, m or p\n");
    exit_status = 1;
    goto END;
  }

#ifdef NAG_COLUMN_MAJOR
  pda = m;
  pdb = p;
  pdc = m;
  pdd = p;
  pdx = n;
#else
  pda = n;
  pdb = n;
  pdc = 1;
  pdd = 1;
  pdx = 1;
#endif

  /* Allocate memory */
  if (!(a = NAG_ALLOC(m * n, Complex)) || !(b = NAG_ALLOC(p * n, Complex)) ||
      !(c = NAG_ALLOC(m, Complex)) || !(d = NAG_ALLOC(p, Complex)) ||
      !(taua = NAG_ALLOC(MIN(m, n), Complex)) ||
      !(taub = NAG_ALLOC(MIN(p, n), Complex)) || !(x = NAG_ALLOC(n, Complex))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read A, B, c and d from data file for the problem
   * min{||c-Ax||_2, x in R^n and B x = d}.
   */
  for (i = 1; i <= m; ++i)
    for (j = 1; j <= n; ++j)
      scanf(" ( %lf , %lf )", &A(i, j).re, &A(i, j).im);
  scanf("%*[^\n]");
  for (i = 1; i <= p; ++i)
    for (j = 1; j <= n; ++j)
      scanf(" ( %lf , %lf )", &B(i, j).re, &B(i, j).im);
  scanf("%*[^\n]");
  for (i = 0; i < m; ++i)
    scanf(" ( %lf , %lf )", &c[i].re, &c[i].im);
  scanf("%*[^\n]");
  for (i = 0; i < p; ++i)
    scanf(" ( %lf , %lf )", &d[i].re, &d[i].im);
  scanf("%*[^\n]");

  /* First compute the generalized RQ factorization of (B,A) as
   * B = (0 R12)*Q,   A = Z*(T11 T12 T13)*Q = T*Q.
   *                        ( 0  T22 T23)
   * where R12, T11 and T22 are upper triangular,
   * using  nag_lapackeig_zggrqf (f08ztc).
   */
  nag_lapackeig_zggrqf(order, p, m, n, b, pdb, taub, a, pda, taua, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_lapackeig_zggrqf (f08ztc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Now, Z^H * (c-Ax) = Z^H * c - T*Q*x, and
   * let f = (f1) = Z^H * (c1) => minimize ||f - T*Q*x||
   *         (f2)         (c2)
   * Compute f using  nag_lapackeig_zunmqr (f08auc), storing result in c
   */
  nag_lapackeig_zunmqr(order, Nag_LeftSide, Nag_ConjTrans, m, 1, MIN(m, n), a,
                       pda, taua, c, pdc, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_lapackeig_zunmqr (f08auc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Putting Q*x = (y1), B * x = d becomes (0 R12) (y1) = d;
   *               (w )                            (w )
   * => R12 * w = d.
   * Solve for w using nag_lapacklin_dtrtrs (f07tec), storing result in d;
   * R12 is (p by p) triangular submatrix starting at B(1,n-p+1).
   */
  nag_lapacklin_ztrtrs(order, Nag_Upper, Nag_NoTrans, Nag_NonUnitDiag, p, 1,
                       &B(1, n - p + 1), pdb, d, pdd, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_lapacklin_ztrtrs (f07tsc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* The problem now reduces to finding the minimum norm of
   * g = (g1) = (f1) - T11*y1 - (T12 T13)*w
   *     (g2)   (f2)          - (T22 T23)*w.
   * Form c1 = f1 - (T12 T13)*w using  nag_blast_zgemv (f16sac).
   */
  alpha = nag_complex_create(-1.0, 0.0);
  beta = nag_complex_create(1.0, 0.0);
  y1rows = n - p;
  nag_blast_zgemv(order, Nag_NoTrans, y1rows, p, alpha, &A(1, n - p + 1), pda,
                  d, 1, beta, c, 1, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_zgemv (f16sac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  /* => now (g1) = c - T11*y1 and ||g1|| = 0 when T11 * y1 = c1.
   * Solve this for y1 using nag_lapacklin_ztrtrs (f07tsc) storing result in c1.
   */
  nag_lapacklin_ztrtrs(order, Nag_Upper, Nag_NoTrans, Nag_NonUnitDiag, y1rows,
                       1, a, pda, c, pdc, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_lapacklin_ztrtrs (f07tsc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* So now Q*x = (y1) is stored in (c1), which is now copied to x.
   *              (w )              (d )
   */
  for (i = 0; i < y1rows; ++i)
    x[i] = nag_complex_create(c[i].re, c[i].im);
  for (i = 0; i < n - y1rows; ++i)
    x[y1rows + i] = nag_complex_create(d[i].re, d[i].im);

  /* Compute x by applying Q inverse using nag_lapackeig_zunmrq (f08cxc). */
  nag_lapackeig_zunmrq(order, Nag_LeftSide, Nag_ConjTrans, n, 1, p, b, pdb,
                       taub, x, pdx, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_lapackeig_zunmrq (f08cxc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* It remains to minimize  ||g2||, g2 = f2 - (T22 T23)*w.
   * Putting w = (y2), gives g2 = f2 - T22*y2 - T23*y3
   *             (y3)
   * [y2 stored in d1, first y2rows of d; y3 stored in d2, next n-m rows of d.]
   *
   * First form T22*y2 using  nag_blast_ztrmv (f16sfc) where y2 is held in d.
   */
  y2rows = MIN(m, n) - y1rows;
  alpha = nag_complex_create(1.0, 0.0);
  nag_blast_ztrmv(order, Nag_Upper, Nag_NoTrans, Nag_NonUnitDiag, y2rows, alpha,
                  &A(n - p + 1, n - p + 1), pda, d, 1, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_ztrmv (f16sfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Then, f2 - T22*y2 (c2 = c2 - d) */
  for (i = 0; i < y2rows; ++i)
    c[y1rows + i] = nag_complex_subtract(c[y1rows + i], d[i]);
  y2rows = m - y1rows;

  if (m < n) {
    y3rows = n - m;
    /* Then g2 = f2 - T22*y2 - T23*y3 (c2 = c2 - T23*d2) */
    alpha = nag_complex_create(-1.0, 0.0);
    nag_blast_zgemv(order, Nag_NoTrans, y2rows, y3rows, alpha,
                    &A(n - p + 1, m + 1), pda, &d[y2rows], 1, beta, &c[y1rows],
                    1, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_blast_zgemv (f16sac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }

  /* Compute ||g|| = ||g2|| = norm(f2 - T22*y2 - T23*y3)
   * using nag_blast_zge_norm (f16uac).
   */
  nag_blast_zge_norm(Nag_ColMajor, Nag_FrobeniusNorm, y2rows, 1, &c[y1rows],
                     y2rows, &rnorm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_blast_zge_norm (f16uac).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Print least squares solution x */
  printf("Constrained least squares solution\n");
  for (i = 0; i < n; ++i)
    printf(" (%7.4f, %7.4f)%s", x[i].re, x[i].im, i % 4 == 3 ? "\n" : "");

  /* Print the square root of the residual sum of squares */
  printf("\nSquare root of the residual sum of squares\n");
  printf("%11.2e\n", rnorm);

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(d);
  NAG_FREE(taua);
  NAG_FREE(taub);
  NAG_FREE(x);

  return exit_status;
}