triuaux.c

matlab中uwb波形优化算法经常会使用的工具包:SeDuMi_1_1R3.

/*
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

*/
#include <string.h>
#include "triuaux.h"
#include "blksdp.h"

/* ************************************************************
   PROCEDURE matperm - Let Y = X(perm,perm)
   INPUT
     x -  n x n input matrix
     perm - length n permutation
     n - order
   OUTPUT
     y - n x n matrix, Y = X(perm,perm)
   ************************************************************ */
void matperm(double *y, const double *x, const int *perm, const int n)
{
  int i,j,inz;
  const double *xj;

  inz = 0;
  for(j = 0; j < n; j++){
    xj = x + perm[j] * n;
    for(i = 0; i < n; i++)
      y[inz++] = xj[perm[i]];
  }
}

/* ************************************************************
   PROCEDURE invmatperm - Let Y(perm,perm) = X
   ************************************************************ */
void invmatperm(double *y, const double *x, const int *perm, const int n)
{
  int i,j,inz;
  double *yj;

  inz = 0;
  for(j = 0; j < n; j++){
    yj = y + perm[j] * n;
    for(i = 0; i < n; i++)
      yj[perm[i]] = x[inz++];
  }
}

/* ************************************************************
   PROCEDURE realltxl: Computes L'*X*L in only n*(n+1)^2/2
   multiplications.
   INPUT
     l - n x n full matrix, tril(l) = L.
     x - n x n full symmetric matrix: instead of X*L, we actually
       compute X'*L.
     n - order.
   OUTPUT
     y - n x n full matrix; tril(Y) = tril(L'*X*L).
   WORK
     xl - length n^2 working vector, to store tril(X'*L).
   ************************************************************ */
void realltxl(double *y, const double *ld, const double *x,
              const int n, double *xld)
{
  int inz, i,j, m, icol, jcol;
/* ------------------------------------------------------------
   Compute xl = tril(X'*L):   n^2 + (n-1)^2 + .. + 1^2 mults.
   For each j, let xj and lj point to x(j,j) and l(j,j); m is
   remaining length of column j, i.e. m=n-j.
   ------------------------------------------------------------ */
  inz = 0; jcol = 0;
  for(j = 0, m = n; j < n; j++, m--, jcol += n+1){
    inz += j;                                 /* go to tril */
    for(i = j, icol = jcol; i < n; i++, icol += n)
      xld[inz++] = realdot(x+icol,ld+jcol,m);   /*xl(i,j) = x(:,i)'*l(j:n,j) */
  }
/* ------------------------------------------------------------
   Compute y = tril(L'*XL):  1*n + 2*(n-1) + ... + n*1 mults.
   ------------------------------------------------------------ */
  inz = 0; jcol = 0;
  for(j = 0; j < n; j++, jcol += n+1){
    inz += j;                                 /* go to tril */
    for(i = j, icol = jcol; i < n; i++, icol += n+1){
      y[inz] = realdot(ld + icol,xld+inz,n-i); /* yij = l(i:n,i)' * xl(:,j) */
      inz++;
    }
  }
}

/* ************************************************************
   PROCEDURE prpiltxl: Computes L'*X*L, with X,L complex, diag(L) real.
   INPUT
     l - n x n full matrix, tril(l) = L. Assume IM diag(L) is all-0.
     x - n x n full Hermitian matrix: instead of X*L, we actually
       compute X'*L.
     n - order.
   OUTPUT
     y - n x n full matrix; tril(Y) = tril(L'*X*L).
   WORK
     xld - length 2 * n^2 working vector, to store tril(X'*L).
   ************************************************************ */
void prpiltxl(double *y,double *ypi, const double *ld,const double *ldpi,
              const double *x,const double *xpi,
              const int n, double *xld)
{
  int inz, i,j, m, icol, jcol;
  double *xldpi;
/* ------------------------------------------------------------
   Partition working array xld into real and imaginary part
   ------------------------------------------------------------ */
  xldpi = xld + SQR(n);
/* ------------------------------------------------------------
   Compute xl = tril(X'*L)        (= tril(XL), since X is Hermitian.)
   For each j, let xj and lj point to x(j,j) and l(j,j); m is
   remaining length of column j, i.e. m=n-j.
   ------------------------------------------------------------ */
  inz = 0; jcol = 0;
  for(j = 0, m = n; j < n; j++, m--, jcol += n+1){
    inz += j;                                 /* go to tril */
    for(i = j, icol = jcol; i < n; i++, icol += n){
/* xl(i,j) = x(:,i)'*l(j:n,j), assume IM l(j,j) == 0 */
      xld[inz] = realdot(x+icol,ld+jcol,m)
        + realdot(xpi+icol+1,ldpi+jcol+1,m-1);
      xldpi[inz] = realdot(x+icol+1,ldpi+jcol+1,m-1)
        - realdot(xpi+icol,ld+jcol,m);
      inz++;
    }
  }
/* ------------------------------------------------------------
   Compute y = tril(L'*XL)
   ------------------------------------------------------------ */
  inz = 0; jcol = 0;
  for(j = 0; j < n; j++, jcol += n+1){
    inz += j;                                 /* go to tril */
    for(i = j, icol = jcol; i < n; i++, icol += n+1){
/* yij = l(i:n,i)' * xl(:,j). Assume IM l(i,i) == 0 */
      y[inz] = realdot(ld + icol,xld+inz,n-i)
        + realdot(ldpi + icol+1,xldpi+inz+1,n-i-1);
      ypi[inz] = realdot(ld + icol,xldpi+inz,n-i)
        - realdot(ldpi+icol+1,xld+inz+1,n-i-1);;
      inz++;
    }
  }
}

/* ************************************************************
   PROCEDURE utmulx - Compute y = triu(U'*XU), U upper triag, in
   1*n + 2*(n-1) + ... + n*1 mults. Only triu(XU) will be used.
   ************************************************************ */
void utmulx(double *y, const double *u, const double *xu, const int n)
{
  int inz,i,j, m;
  const double *ui;

  inz = 0;
  for(j = 0, m = n; j < n; j++, xu += n){
    for(i = 0, ui = u; i <= j; i++, ui += n)
      y[inz++] = realdot(ui,xu,i+1);     /* yij = u(0:i,i)' * xu(0:i,j) */
    inz += --m;                          /* skip tril */
  }
}

/* ************************************************************
   PROCEDURE prpiutmulx - Compute y = triu(U'*XU), U upper triag, in
   4*(1*n + 2*(n-1) + ... + n*1) - n(n-1) mults. Only triu(XU) will be used.
   CAUTION: assumes that diag(imag(U)) == 0.
   ************************************************************ */
void prpiutmulx(double *y, double *ypi, const double *u, const double *upi,
                const double *xu, const double *xupi, const int n)
{
  int inz,i,j, m, icol;

  inz = 0;
  for(j = 0, m = n; j < n; j++, xu += n, xupi += n){
    for(i = 0, icol = 0; i <= j; i++, icol += n){
/* yij = u(0:i,i)' * xu(0:i,j) where IM(u(i,i))==0*/
      y[inz] = realdot(u+icol,xu,i+1) + realdot(upi + icol, xupi,i);
      ypi[inz] = realdot(u+icol,xupi,i+1) - realdot(upi + icol, xu, i);
      inz++;
    }
    inz += --m;                          /* skip tril */
  }
}

/* ************************************************************
   PROCEDURE realutxu: Computes U'*X*U in only n*(n+1)^2/2
   multiplications.
   INPUT
     u - n x n full matrix, triu(u) = U.
     x - n x n full symmetric matrix: instead of X*U, we actually
       compute X'*U.
     n - order.
   OUTPUT
     y - n x n full matrix; triu(Y) = triu(U'*X*U).
   WORK
     xu - length n^2 working vector, to store triu(X'*U).
   ************************************************************ */
void realutxu(double *y, const double *u, const double *x,
              const int n, double *xu)
{
  int inz, i,j, m, icol,jcol;
/* ------------------------------------------------------------
   Compute xu = triu(X'*U):   1^2 + 2^2 + ... + (n-1)^2 + n^2 mults.
   For each j, let uj point to u(0,j); m is
   remaining length of column j, i.e. m=j.
   ------------------------------------------------------------ */
  inz = 0; jcol = 0;
  for(j = 0, m = n; j < n; j++, jcol += n){
    for(i = 0, icol = 0; i <= j; i++, icol += n)
      xu[inz++] = realdot(x+icol,u+jcol,j+1);  /*xu(i,j) = x(:,i)'*u(0:j,j) */
    inz += --m;                                /* skip tril */
  }
/* ------------------------------------------------------------
   Compute y = triu(U'*XU):  1*n + 2*(n-1) + ... + n*1 mults.
   ------------------------------------------------------------ */
  utmulx(y, u, xu, n);
}

/* ************************************************************
   PROCEDURE prpiutxu: Computes U'*X*U with complex data.
   INPUT
     u - n x n full matrix, triu(u) = U.
     x - n x n full symmetric matrix: instead of X*U, we actually
       compute X'*U.
     n - order.
   OUTPUT
     y - n x n full matrix; triu(Y) = triu(U'*X*U).
   WORK
     xu - length n^2 working vector, to store triu(X'*U).
   ************************************************************ */
void prpiutxu(double *y, double *ypi, const double *u, const double *upi,
              const double *x, const double *xpi,
              const int n, double *xu)
{
  int inz, i,j, m, icol,jcol;
  double *xupi;
/* ------------------------------------------------------------
   Partition working array xu into real and imaginary part
   ------------------------------------------------------------ */
  xupi = xu + SQR(n);
/* ------------------------------------------------------------
   Compute xu = triu(X'*U)   ( =triu(XU), since X is Hermitian.)
   For each j, let uj point to u(0,j); m is
   remaining length of column j, i.e. m=j.
   ------------------------------------------------------------ */
  inz = 0; jcol = 0;
  for(j = 0, m = n; j < n; j++, jcol += n){
    for(i = 0, icol = 0; i <= j; i++, icol += n){
/*xu(i,j) = x(:,i)'*u(0:j,j), assume IM u(j,j) == 0 */
      xu[inz] = realdot(x+icol,u+jcol,j+1) + realdot(xpi+icol,upi+jcol,j);
      xupi[inz] = realdot(x+icol,upi+jcol,j) - realdot(xpi+icol,u+jcol,j+1);
      inz++;
    }
    inz += --m;                                /* skip tril */
  }
/* ------------------------------------------------------------
   Compute y = triu(U'*XU)
   ------------------------------------------------------------ */
  prpiutmulx(y, ypi, u, upi, xu, xupi, n);
}

/* ************************************************************
   PROCEDURE partfwsolve -- Solve y from U(1:n,1:n)'*y(1:n) = x(1:n)
      where U is upper-triangular, and of order m >= n.
   INPUT
     u - m x m full matrix; only its upper-triangular entries get used
     x - length n vector.
     m - order of u, m >= n
     n - order of x
   OUTPUT
     y - length n vector, y = (U'\x)_{1:n}.
   ************************************************************ */
void partfwsolve(double *y, const double *u, const double *x,
                const int m, const int n)
{
  int k;
/* ------------------------------------------------------------
   The first equation, u(:,1)'*y=x(1), yields y(1) = x(1)/u(1,1).
   For k = 2:n, we solve
   u(1:k-1,k)'*y(1:k-1) + u(k,k)*y(k) = x(k).
   ------------------------------------------------------------ */
  y[0] = x[0] / u[0];
  u += m;                             /* done with the first column of u */
  for(k = 1; k < n; k++, u += m)
    y[k] = (x[k] - realdot(y,u,k)) / u[k];
}

/* Complex case. Assume IM diag(u) is all-0.
 Solves U' * y = x */
void prpipartfwsolve(double *y,double *ypi, const double *u,const double *upi,
                     const double *x,const double *xpi,
                     const int m, const int n)
{
  int k;
/* ------------------------------------------------------------
   The first equation, u(:,1)'*y=x(1), yields y(1) = x(1)/u(1,1).
   For k = 2:n, we solve
   u(1:k-1,k)'*y(1:k-1) + u(k,k)*y(k) = x(k).
   Assume that IM u(k,k) == 0.
   ------------------------------------------------------------ */
  y[0] = x[0] / u[0];
  ypi[0] = xpi[0] / u[0];
  u += m;                             /* done with the first column of u */
  upi += m;
  for(k = 1; k < n; k++, u += m, upi+=m){
/* y(k) = {x(k) - u(1:k-1,k)'*y(1:k-1)} / ukk */
    y[k] = (x[k] - realdot(u,y,k) - realdot(upi,ypi,k)) / u[k];
    ypi[k] = (xpi[k] - realdot(u,ypi,k) + realdot(upi,y,k)) / u[k];
  }
}

/* ************************************************************
   PROCEDURE partbwsolve -- Solve y from L(j:m,j:m)'*y(j:m) = x(j:m)
      where L is lower-triangular, and of order m > j.
   INPUT
     l - m x m full matrix; only its lower-triangular entries get used
     x - length n vector.
     m - order of l, m > j
     j - last entry of x to be computed (going backwards). If j = 0, then
       all m entries are computed, otherwise only the m-j bottom ones.
       0 <= j < m.
   OUTPUT
     y - length m vector, y[j:m-1] = (L'\x)_{[j:m-1]}.
   ************************************************************ */
void partbwsolve(double *y, const double *l, const double *x,
                 const int m, const int j)
{
  int k;
  const double *lk;
  double ldoty;
/* ------------------------------------------------------------
   The last equation, l(:,m)'*y=x(m), yields y(m) = x(m)/l(m,m).
   For k = m-1:j, we solve
   l(k+1:m,k)'*y(k+1:m) + l(k,k)*y(k) = x(k).
   ------------------------------------------------------------ */
  k = m-1;
  lk = l + m*k;
  ldoty = 0.0;
  for(k = m-1; k > j; k--){
    y[k] = (x[k] - ldoty) / lk[k];
    lk -= m;
    ldoty = realdot(y+k,lk+k,m-k);       /* ldoty = l(k:m,k-1)'*y(k:m) */
  }
  y[j] = (x[j] - ldoty) / lk[j];
}