function [cmps, cmpa, rc] = mfs_statrespc(cmp, opts)

# usage: [cmps, cmpa, rc] = mfs_statrespc(cmp, opts)
#
# Input  cmp     Structure with data of aeroelastic component
#        opts    Structure with options
#
# Output cmps    Structure with data of solid component
#        cmpa    Structure with data of aerodynamic component
#        rc      Return code: 0 = no errors
#                             1 = errors
#
# The function computes the static aeroelastic response of a
# structure without rigid body modes. The results are stored in the
# solid and aerodynamic components respectively.
#
# --------------------------------------------------------------------

  rc   = 0;
  cmps = struct();
  cmpa = struct();

# Check arguments

  if (nargin != 2 || nargout != 3)
     print_usage();
  end

# Check component status

  if (! isfield(cmp.solid, "stiff"))
     printf("*E* mfs_statresp: stiffness matrix missing\n");
     rc = 1;
  end
  if (! cmp.aero.config.nconf)
     printf("*E* mfs_statresp: no configurations defined\n");
     rc = 1;
  end
  if (! isfield(cmp, "splines"))
     printf("*E* mfs_statresp: definition of splines missing\n");
     rc = 1;
  end
  if (! isfield(cmp.splines, "Shg"))
     printf("*E* mfs_statresp: spline matrices missing\n");
     rc = 1;
  end

  if (rc) return; end

# Initialize output

  cmps = cmp.solid;
  cmpa = cmp.aero;

# Downwash due to camber and flaps

  dw      = mfs_vlmdw(cmp.aero.panels, cmp.aero.controls,
                      cmp.aero.config);  
  confmat = mfs_confmat(cmp.aero.config);
  wr      = dw * confmat;

# Matrix C to compute downwash from vortices

  symy = isfield(cmp.aero, "symy");
  if (symy)
     ys = cmp.aero.symy;
  else
     ys = 0;
  end

  C = mfs_vlmmat(cmp.aero.ls, cmp.aero.panels, symy, ys);

# Combine normal wash with normal wash due to elastic motion

  w   = [wr, cmp.splines.D1h];
  lw  = columns(w);
  lwr = columns(wr);
  lw1 = lwr + 1;

# Vortex matrix

  Vxmat = C \ w;

# Matrix G to compute forces from vortices

  G = 2 * diag(cmp.aero.panels.width);

# Spline matrix to local displacements

  Shl = mfs_matpartc(cmp.splines.Shg, cmp.solid.dofs);

# Aerodynamic matrices

  Svh = G * cmp.splines.Svh;
  Vxh = Svh' * Vxmat;
  qrg = cmp.splines.Shg' * Vxh(:, 1 : lwr);
  qrl = mfs_matpartr(qrg, cmp.solid.dofs);
  Qhh = Vxh(:, lw1 : lw);

# Stiffness matrix

  Kll = mfs_matpart(cmp.solid.stiff.K, cmp.solid.dofs, [1, 0]);

# Loads on solid

  % to be inserted here

# Initialize displacements

  nconf = cmp.aero.config.nconf;
  disp  = zeros(cmp.solid.dofs.ndofg, nconf, "double");

# Compute common matrices

  uQ  = Kll \ [qrl, Shl' * Qhh];
  Qlh = uQ(:, nconf + 1 : end);

  dofl = cmp.solid.dofs.dofl;

# Branch according to method

  switch opts.method

  case "direct"
     Ihh = eye(rows(Qhh));
     for n = 1 : nconf
         qdyn = cmp.aero.config.qdyn(n);
         Rhh  = Ihh - (qdyn * Shl) * Qlh;
         u0   = qdyn * uQ(:, n);
         u0  += Qlh * (Rhh \ (Shl * (qdyn * u0)));
         disp(dofl, n) = u0;
     end 

  case "bicg"
     Afn = @(x, type, qdyn) Amat(x, type, qdyn, Qlh, Shl);
     idn = @(x, y, z) x;  % avoid bug in 7.1.0
     for n = 1 : nconf
         niter = 0;
         qdyn  = cmp.aero.config.qdyn(n);
         u0    = qdyn * uQ(:, n);
         [disp(dofl, n), flag, ~, niter] = ...
              bicg(Afn, u0, opts.tol, opts.maxit, idn, idn,
                   u0, qdyn);
         printf("           Configuration %2.0d: ", n);
         if (! flag)
            printf("solution converged within %2.0d iterations\n",
                   niter);
         else
            printf("solution failed, bicg flag = %2.0d, iter = %3.0d\n", 
                   flag, niter);
            rc = 1;
         end
     end

  otherwise
     printf("*E* mfs_statresp: unknown method \"%s\"\n",
            opts.method);
     rc = 1;

  end

  if (rc) return; end

  if ((ndofd = cmp.solid.dofs.ndofd))
     dofd = cmp.solid.dofs.dofd;
     disp(dofd, :) = cmp.solid.dofs.C(dofd, :) * disp;
  end

# Vortices

  dh    = Shl * disp(dofl, :);
  Gamma = Vxmat(:, 1 : lwr) + Vxmat(:, lw1 : lw) * dh;

# Compute reaction loads and internal constraint loads

  dofp     = cmp.solid.dofs.dofp;
  mxdofpnt = cmp.solid.dofs.mxdofpnt;

  qdyn = diag([cmp.aero.config.qdyn]);
  h    = cmp.splines.Shg * disp;

  FE = cmp.solid.stiff.K * disp;
  FA = qrg + cmp.splines.Shg' * (Qhh * h);
  F  = FE - FA * qdyn;

  if (ndofd)
     Fp = cmp.solid.dofs.C(:, dofp)' * F;
     F(dofp, :) -= Fp;
  else
     Fp = F(dofp, :);
  end

  reac = mfs_mat2nodedata(Fp, dofp, mxdofpnt, 4);
  reac = mfs_resnodedata(reac, cmp.solid.nodes);

  if (ndofd)
     dofda = unique([cmp.solid.dofs.dofa; dofd]);
     icsl  = mfs_mat2nodedata(F(dofda, :), dofda, mxdofpnt, 4);
     icsl  = mfs_resnodedata(icsl, cmp.solid.nodes);
  end

# Store results

  resps = struct("nofldc", nconf, "disp", disp, "reac", reac);

  if (ndofd)
     resps.icsl = icsl;
  end

  respa = struct("nconf", nconf, "G", Gamma);
  respa.disp = cmp.splines.Snh * (cmp.splines.Shg * disp);

  cmps.statresp = resps;
  cmpa.statresp = respa;

end

function y = Amat(x, type, qdyn, Qlh, Shl)

# This function computes the matrix product 
#       (Ill - qdyn * Kll \ Qll) * x.
# It is used by the iterative solver.
#
# ---------------------------------------------------------------------

  switch type
  case "notransp"
     y = Qlh * (Shl * x);
  case "transp"
     y = Shl' * (Qlh' * x);
  end

  y = x - qdyn * y;

end
