関連情報
ホーム > 製品 > NAG数値計算ライブラリ > サンプルソースコード集 > 2次計画問題(密な) (C言語/C++)

2次計画問題(密な)

C言語によるサンプルソースコード
使用関数名:nag_opt_qp (e04nfc)

Keyword: 2次計画問題, 密な

概要

本サンプルは2次計画問題(密な)を解くC言語によるサンプルプログラムです。 本サンプルは以下に示される制約条件を満たす2次関数を最小化する解を求めて出力します。

2次計画問題のデータ 

※本サンプルはNAG Cライブラリに含まれる関数 nag_opt_qp() のExampleコードです。本サンプル及び関数の詳細情報は nag_opt_qp のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

オプション・パラメータ

(本関数の詳細はnag_opt_qp のマニュアルページを参照)

このデータをダウンロード
nag_opt_qp (e04nfc) Example Program Optional Parameters

Following options for e04nfc are read by e04xyc.

begin e04nfc

 fmax_iter = 30 /* Set maximum number of iterations in feasiblity phase */
 max_iter = 50  /* Set maximum total number of iterations */

end

  • 1行目はタイトル行で読み飛ばされます。
  • 5〜10行目にe04xycによって読み込まれる、e04nfcに渡すオプショナル引数を指定しています。
    fmax_iter 実行可能フェーズでの最大反復数。
    max_iter 最大反復数。

入力データ

(本関数の詳細はnag_opt_qp のマニュアルページを参照)

このデータをダウンロード
nag_opt_qp (e04nfc) Example Program Data
Linear term of f(x), c.
-0.02 -0.2 -0.2 -0.2 -0.2 0.04 0.04
Linear constraint matrix, A.
1.0  1.0  1.0  1.0  1.0  1.0  1.0
0.15 0.04 0.02 0.04 0.02 0.01 0.03
0.03 0.05 0.08 0.02 0.06 0.01 0.0
0.02 0.04 0.01 0.02 0.02 0.0  0.0
0.02 0.03 0.0  0.0  0.01 0.0  0.0
0.70 0.75 0.80 0.75 0.80 0.97 0.0
0.02 0.06 0.08 0.12 0.02 0.01 0.97
Lower bounds
-0.01 -0.1 -0.01 -0.04 -0.1 -0.01 -0.01
-0.13 -1.0e21 -1.0e21 -1.0e21 -1.0e21 -0.0992 -0.003
Upper bounds
 0.01  0.15 0.03  0.02  0.05 1.0e21 1.0e21
-0.13 -0.0049 -0.0064 -0.0037 -0.0012  1.0e21  0.002
Initial estimate of x
-0.01 -0.03 0.0 -0.01 -0.1 0.02 0.01

  • 1行目はタイトル行で読み飛ばされます。
  • 3行目にf(x)の線形項(cvec)を指定しています。
  • 5〜11行目には線形制約行列Aの要素を指定しています。
  • 13〜14行目には変数と制約の下限(bl)を指定しています。
  • 16〜17行目には変数と制約の上限(bu)を指定しています。
  • 19行目にはxの初期推定値を指定しています。

出力結果

(本関数の詳細はnag_opt_qp のマニュアルページを参照)

この出力例をダウンロード
nag_opt_qp (e04nfc) Example Program Results

Optional parameter setting for e04nfc.
--------------------------------------

Option file: e04nfce.opt

fmax_iter set to 30
max_iter set to 50

Parameters to e04nfc
--------------------

Linear constraints............   7    Number of variables...........   7

prob....................   Nag_QP2    start...................  Nag_Cold
ftol.................... 1.05e-008    reset_ftol..............         5
rank_tol................ 1.11e-014    crash_tol............... 1.00e-002
fcheck..................        50    max_df..................         7
inf_bound............... 1.00e+021    inf_step................ 1.00e+021
fmax_iter...............        30    max_iter................        50
hrows...................         7    machine precision....... 1.11e-016
optim_tol............... 1.72e-013    min_infeas.............. Nag_FALSE
print_level......... Nag_Soln_Iter
outfile.................    stdout

Memory allocation:
state...................       Nag
ax......................       Nag    lambda..................       Nag

Results from e04nfc:
-------------------

  Itn Jdel  Jadd   Step    Ninf  Sinf/Obj    Bnd  Lin  Nart  Nrz  Norm Gz

   0   0     0   0.0e+000   3  1.0380e-001    3    4    0     0  0.00e+000
   1   9 U  13 L 4.1e-002   1  3.0000e-002    3    4    0     0  0.00e+000
   2  12 U   4 L 4.2e-002   0  0.0000e+000    4    3    0     0  0.00e+000

 Itn 2  -- Feasible point found.
   2   0     0   0.0e+000   0  4.5800e-002    4    3    0     0  0.00e+000
   3   5 L  14 L 1.3e-001   0  4.1616e-002    3    4    0     0  0.00e+000
   4  11 U   0   1.0e+000   0  3.9362e-002    3    3    0     1  2.78e-017
   5   3 L  10 U 4.1e-001   0  3.7589e-002    2    4    0     1  1.19e-002
   6   0     0   1.0e+000   0  3.7554e-002    2    4    0     1  1.04e-017
   7   4 L   0   1.0e+000   0  3.7032e-002    1    4    0     2  5.56e-017

Final solution:

  Varbl State    Value      Lower Bound  Upper Bound    Lagr Mult    Residual

 V   1    LL -1.00000e-002 -1.0000e-002  1.0000e-002   4.700e-001  0.000e+000
 V   2    FR -6.98646e-002 -1.0000e-001  1.5000e-001   0.000e+000  3.014e-002
 V   3    FR  1.82592e-002 -1.0000e-002  3.0000e-002   0.000e+000  1.174e-002
 V   4    FR -2.42608e-002 -4.0000e-002  2.0000e-002   0.000e+000  1.574e-002
 V   5    FR -6.20056e-002 -1.0000e-001  5.0000e-002   0.000e+000  3.799e-002
 V   6    FR  1.38054e-002 -1.0000e-002   None      0.000e+000  2.381e-002
 V   7    FR  4.06650e-003 -1.0000e-002   None      0.000e+000  1.407e-002

   LCon  State    Value     Lower Bound  Upper Bound    Lagr Mult    Residual

 L   1    EQ -1.30000e-001 -1.3000e-001 -1.3000e-001  -1.908e+000  2.776e-017
 L   2    FR -5.87990e-003     None     -4.9000e-003   0.000e+000  9.799e-004
 L   3    UL -6.40000e-003     None     -6.4000e-003  -3.144e-001  0.000e+000
 L   4    FR -4.53732e-003     None     -3.7000e-003   0.000e+000  8.373e-004
 L   5    FR -2.91600e-003     None     -1.2000e-003   0.000e+000  1.716e-003
 L   6    LL -9.92000e-002 -9.9200e-002   None      1.955e+000  0.000e+000
 L   7    LL -3.00000e-003 -3.0000e-003  2.0000e-003   1.972e+000 -4.337e-019

Exit after 7 iterations.

Optimal QP solution found.

Final QP objective value =  3.7031646e-002


A run of the same example with a warm start:

Parameters to e04nfc
--------------------

Linear constraints............   7    Number of variables...........   7

prob....................   Nag_QP2    start...................  Nag_Warm
ftol.................... 1.05e-008    reset_ftol..............         5
rank_tol................ 1.11e-014    crash_tol............... 1.00e-002
fcheck..................        50    max_df..................         7
inf_bound............... 1.00e+021    inf_step................ 1.00e+021
fmax_iter...............        30    max_iter................        50
hrows...................         7    machine precision....... 1.11e-016
optim_tol............... 1.72e-013    min_infeas.............. Nag_FALSE
print_level.........      Nag_Soln
outfile.................    stdout

Memory allocation:
state...................       Nag
ax......................       Nag    lambda..................       Nag

Final solution:

  Varbl State    Value      Lower Bound  Upper Bound    Lagr Mult    Residual

 V   1    LL -1.00000e-002 -1.0000e-002  1.0000e-002   4.700e-001  0.000e+000
 V   2    FR -6.98646e-002 -1.0000e-001  1.5000e-001   0.000e+000  3.014e-002
 V   3    FR  1.82592e-002 -1.0000e-002  3.0000e-002   0.000e+000  1.174e-002
 V   4    FR -2.42608e-002 -4.0000e-002  2.0000e-002   0.000e+000  1.574e-002
 V   5    FR -6.20056e-002 -1.0000e-001  5.0000e-002   0.000e+000  3.799e-002
 V   6    FR  1.38054e-002 -1.0000e-002   None      0.000e+000  2.381e-002
 V   7    FR  4.06650e-003 -1.0000e-002   None      0.000e+000  1.407e-002

   LCon  State    Value     Lower Bound  Upper Bound    Lagr Mult    Residual

 L   1    EQ -1.30000e-001 -1.3000e-001 -1.3000e-001  -1.908e+000  0.000e+000
 L   2    FR -5.87990e-003     None     -4.9000e-003   0.000e+000  9.799e-004
 L   3    UL -6.40000e-003     None     -6.4000e-003  -3.144e-001  0.000e+000
 L   4    FR -4.53732e-003     None     -3.7000e-003   0.000e+000  8.373e-004
 L   5    FR -2.91600e-003     None     -1.2000e-003   0.000e+000  1.716e-003
 L   6    LL -9.92000e-002 -9.9200e-002   None      1.955e+000 -1.388e-017
 L   7    LL -3.00000e-003 -3.0000e-003  2.0000e-003   1.972e+000 -2.602e-018

Exit after 0 iterations.

Optimal QP solution found.

Final QP objective value =  3.7031646e-002


  • 6行目にオプション・ファイルe04nfce.optが読み込まれたことが示されています。
  • 8行目にfmax_iter が30に設定されたことが示されています。
  • 9行目にmax_iter が30に設定されたことが示されています。
  • 14行目に線形制約の数、変数の数が出力されています。
  • 16〜25行目に以下に示すプログラムのオプション引数が出力されています。
    prob 最小化される目的関数の種類。"Nag_QP2"は2次計画問題であることを意味します。
    start 初期のワーキングセットがどのように選ばれるかを示しています。"Nag_Cold"は変数や制約の値に基づいて初期のワーキングセットが選ばれていることを意味します。
    ftol 実行可能解での許容可能な制約違反の最大値。
    reset_ftol 問題を解くためにこの反復数を超えた反復が必要な場合に、新たに反復が開始されます。
    fcheck この反復回数の時に現在の解がワーキングセットにある制約を満たしているか確認されます。
    max_df 縮小へシアン(Reduced Hessian)の三角因子Rに割り当てられた記憶域の上限。
    crash_tol 初期のワーキングセットに含まれる不等式制約が下限/上限のこの範囲内にあることを示しています。
    inf_bound 無限の境界。この値以上の上限は+∞と見なされます。またこの値を負に反転させた値以下の下限は−∞と見なされます。
    inf_step 無限解へのステップと見なされる変数の変化の大きさ。
    fmax_iter 実行可能フェーズでの最大反復数。
    max_iter 最適化フェーズの最大反復数。
    hrows QP目的関数の二次項Hの行数。
    machine precision マシンの精度。
    optim_tol 境界値や制約に正しい符号がついているかを判断するのに使用される許容値。
    min_infeas 制約に対して実行可能解がない場合、実行不可能性を最小化すべきかを指定するフラグ。 "Nag_FALSE" を指定した場合、問題が実行不可能とわかり次第ルーチンは終了します。
    print_level 結果出力のレベル。"Nag_Soln_Iter"は最終解と各反復の結果を出力することを意味します。
    outfile 結果が出力されるファイル名。"stdout"は標準出力を意味します。
  • 27〜29行目には、オプション引数state、ax、lambdaに必要なメモリがNAGライブラリの関数によって自動的に割り当てられたことを示しています。
  • 34〜46行目には以下が出力されています。また40行目に2回目の反復で実行可能点が見つかったことが示されています。
    Itn 反復回数。
    Jdel ワーキングセットから削除される制約のインデックス。"U"は上限、"L"は下限を意味します。
    Jadd ワーキングセットに追加される制約のインデックス。"U"は上限、"L"は下限を意味します。
    Step 探索方向に進んだステップ幅。
    Ninf 違反した制約の数。
    Sinf/Obj 制約違反の大きさの加重和(xが実行可能でない場合)もしくは目的関数の値(xが実行可能な場合)。
    Bnd 現在のワーキングセットの簡易境界制約の数。
    Lin 現在のワーキングセットの一般線形制約の数。
    Nart ワーキングセットの人為的制約の数。
    Nrz 目的関数が最小化される部分空間の次元数。
    Norm Gz 縮小勾配のユークリッド・ノルム。
  • 50〜58行目には以下が出力されています。
    Varbl 変数を示す"V"、インデックス。
    State 変数の状態。FRは下限も上限もワーキングセットにはないことを意味し、LLは下限であることを意味します。
    Value 最後の反復での変数の値。
    Lower Bound 下限。
    Upper Bound 上限。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
    Residual 残差(変数の値と上限/下限の近いほうとの差)。
  • 60〜68行目には以下が出力されています。
    LCon 線形制約を示す"L"、インデックス。
    State 制約の状態。EQは固定制約を意味し、LLは下限であることを意味し、ULは上限であることを意味し、FRは下限も上限もワーキングセットにはないことを意味します。
    Value 最後の反復での制約の値。
    Lower Bound 下限。
    Upper Bound 上限。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
    Residual 残差(制約の値と上限/下限の近いほうとの差)。
  • 70行目に7回の反復を行ったことが示されています。
  • 72行目にQPの最適解が見つかったことが示されています。
  • 74行目にQPの最終的な目的関数値が出力されています。
  • 77行目にwarm startで同じQP問題を解くことが示されています。
  • 82行目に線形制約の数、変数の数が出力されています。
  • 84〜93行目にプログラムのオプション引数が出力されています。
  • 95〜97行目には、オプション引数state、ax、lambdaに必要なメモリがNAGライブラリの関数によって自動的に割り当てられたことを示しています。
  • 101〜109行目には以下が出力されています。
    Varbl 変数を示す"V"、インデックス。
    State 変数の状態。FRは下限も上限もワーキングセットにはないことを意味し、LLは下限であることを意味します。
    Value 最後の反復での変数の値。
    Lower Bound 下限。
    Upper Bound 上限。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
    Residual 残差(変数の値と上限/下限の近いほうとの差)。
  • 111〜119行目には以下が出力されています。
    LCon 線形制約を示す"L"、インデックス。
    State 制約の状態。EQは固定制約を意味し、LLは下限であることを意味し、ULは上限であることを意味し、FRは下限も上限もワーキングセットにはないことを意味します。
    Value 最後の反復での制約の値。
    Lower Bound 下限。
    Upper Bound 上限。
    Lagr Mult 関連する境界制約のラグランジュ乗数。
    Residual 残差(制約の値と上限/下限の近いほうとの差。)
  • 121行目に反復せずに終了したことが示されています。
  • 123行目にQPの最適解が見つかったことが示されています。
  • 125行目にQPの最終的な目的関数値が出力されています。

ソースコード

(本関数の詳細はnag_opt_qp のマニュアルページを参照)

※本サンプルソースコードはNAG数値計算ライブラリ(Windows, Linux, MAC等に対応)の関数を呼び出します。
サンプルのコンパイル及び実行方法


このソースコードをダウンロード
/* nag_opt_qp (e04nfc) Example Program.
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 * Mark 6 revised, 2000.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

#include <nag_example_file_io.h>
#include <nag.h>
#include <stdio.h>
#include <string.h>
#include <nag_stdlib.h>
#include <nag_string.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL qphess1(Integer n, Integer jthcol, double h[],
                             Integer tdh, double x[], double hx[],
                             Nag_Comm *comm);
static void NAG_CALL qphess2(Integer n, Integer jthcol, double h[],
                             Integer tdh, double x[], double hx[],
                             Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

#define A(I, J) a[(I) *tda + J]
#define H(I, J) h[(I) *tdh + J]

int main(int argc, char *argv[])
{
  FILE        *fpin, *fpout;
  char        *optionsfile = 0;
  char        *outfile = 0;
  Nag_Boolean print;
  Integer     exit_status = 0, i, j, n, nbnd, nclin, tda, tdh;
  Nag_E04_Opt options;
  double      *a = 0, *bl = 0, *bu = 0, *cvec = 0, *h = 0, objf, *x = 0;
  NagError    fail;

  INIT_FAIL(fail);


  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  (void) nag_example_file_io(argc, argv, "-options", &optionsfile);
  (void) nag_example_file_io(argc, argv, "-nag_write", &outfile);
  if (!outfile)
    {
      outfile = NAG_ALLOC(7, char);
      strncpy(outfile, "stdout", strlen("stdout")+1);
    }

  fprintf(fpout, "nag_opt_qp (e04nfc) Example Program Results\n");

  fscanf(fpin, " %*[^\n]"); /* Skip heading in data file */

  /* Set the actual problem dimensions.
   * n = the number of variables.
   * nclin = the number of general linear constraints (may be 0).
   */
  n = 7;
  nclin = 7;
  if (n > 0 && nclin >= 0)
    {
      nbnd = n + nclin;
      if (!(x = NAG_ALLOC(n, double)) ||
          !(cvec = NAG_ALLOC(n, double)) ||
          !(a = NAG_ALLOC(nclin*n, double)) ||
          !(h = NAG_ALLOC(n*n, double)) ||
          !(bl = NAG_ALLOC(nbnd, double)) ||
          !(bu = NAG_ALLOC(nbnd, double)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
      tdh = n;
    }
  else
    {
      fprintf(fpout, "Invalid n or nclin.\n");
      exit_status = 1;
      return exit_status;
    }
  /* cvec   = the coefficients of the explicit linear term of f(x).
   * a      = the linear constraint matrix.
   * bl     = the lower bounds on x and A*x.
   * bu     = the upper bounds on x and A*x.
   * x      = the initial estimate of the solution.
   */

  /* Read the coefficients of the explicit linear term of f(x). */
  fscanf(fpin, " %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < n; ++i)
    fscanf(fpin, "%lf", &cvec[i]);

  /* Read the linear constraint matrix A. */
  fscanf(fpin, " %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < nclin; ++i)
    for (j = 0; j < n; ++j)
      fscanf(fpin, "%lf", &A(i, j));

  /* Read the bounds. */
  nbnd = n + nclin;
  fscanf(fpin, " %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < nbnd; ++i)
    fscanf(fpin, "%lf", &bl[i]);
  fscanf(fpin, " %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < nbnd; ++i)
    fscanf(fpin, "%lf", &bu[i]);

  /* Read the initial estimate of x. */
  fscanf(fpin, " %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < n; ++i)
    fscanf(fpin, "%lf", &x[i]);

  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options); /* Initialise options structure */
  strncpy(options.outfile, outfile, strlen(outfile)+1);
  /* Set one option directly
   * Bounds  >=   inf_bound will be treated as plus  infinity.
   * Bounds  <=  -inf_bound will be treated as minus infinity.
   */
  options.inf_bound = 1.0e21;

  /* Read remaining option values from file */
  print = Nag_TRUE;
  /* nag_opt_read (e04xyc).
   * Read options from a text file
   */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_read("e04nfc", optionsfile, &options, print, options.outfile, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_opt_read (e04xyc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Solve the problem from a cold start.
   * The Hessian is defined implicitly by function qphess1.
   */
  /* nag_opt_qp (e04nfc), see above. */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, (double *) 0, tdh,
             qphess1, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_opt_qp (e04nfc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* The following is for illustrative purposes only. We do a warm
   * start with the final working set of the previous run.
   * This time we store the Hessian explicitly in h[][], and use
   * the corresponding function qphess2().
   * Only the final solution from the results is printed.
   */
  fprintf(fpout, "\nA run of the same example with a warm start:\n");

  options.start = Nag_Warm;
  options.print_level = Nag_Soln;

  for (i = 0; i < n; ++i)
    {
      for (j = 0; j < n; ++j) H(i, j) = 0.0;
      if (i <= 4) H(i, i) = 2.0;
      else H(i, i) = -2.0;
    }
  H(2, 3) = 2.0;
  H(3, 2) = 2.0;
  H(5, 6) = -2.0;
  H(6, 5) = -2.0;

  /* Solve the problem again. */
  /* nag_opt_qp (e04nfc), see above. */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, h, tdh,
             qphess2, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_opt_qp (e04nfc).\n%s\n", fail.message);
      exit_status = 1;
    }
  /* Free memory allocated by nag_opt_qp (e04nfc) to pointers in options */
  /* nag_opt_free (e04xzc).
   * Memory freeing function for use with option setting
   */
  nag_opt_free(&options, "all", &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (x) NAG_FREE(x);
  if (cvec) NAG_FREE(cvec);
  if (a) NAG_FREE(a);
  if (h) NAG_FREE(h);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (optionsfile) NAG_FREE(optionsfile);
  if (outfile) NAG_FREE(outfile);

  return exit_status;
}

static void NAG_CALL qphess1(Integer n, Integer jthcol, double h[],
                             Integer tdh, double x[], double hx[],
                             Nag_Comm *comm)
{
  /* In this version of qphess the Hessian matrix is implicit.
   * The array h[] is not accessed. There is no special coding
   * for the case jthcol > 0.
   */

  hx[0] = 2.0*x[0];
  hx[1] = 2.0*x[1];
  hx[2] = 2.0*(x[2] + x[3]);
  hx[3] = hx[2];
  hx[4] = 2.0*x[4];
  hx[5] = -2.0*(x[5] + x[6]);
  hx[6] = hx[5];
} /* qphess1 */

#undef H

static void NAG_CALL qphess2(Integer n, Integer jthcol, double h[],
                             Integer tdh, double x[], double hx[],
                             Nag_Comm *comm)
{
  /* In this version of qphess, the matrix H is stored in h[]
   * as a full two-dimensional array.
   */

#define H(I, J) h[(I) *tdh + (J)]

  Integer i, j;

  if (jthcol != 0)
    {
      /* Special case -- extract one column of  H. */
      j = jthcol - 1;
      for (i = 0; i < n; ++i)
        hx[i] = H(i, j);
    }
  else
    {
      /* Normal Case. */
      for (i = 0; i < n; ++i) hx[i] = 0.0;

      for (i = 0; i < n; ++i)
        for (j = 0; j < n; ++j)
          hx[i] += H(i, j)*x[j];
    }
} /* qphess2 */


Results matter. Trust NAG.

Privacy Policy | Trademarks