/* nag_real_partial_svd (f02wgc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 9, 2009.
*/
/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagf02.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL av(Integer *iflag, Integer m, Integer n, const double x[],
double ax[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
int main(void)
{
/*Integer scalar and array declarations */
Integer       exit_status = 0;
Integer       i, m, n, nconv, ncv, nev;
Integer       pdu, pdv;
Nag_Comm      comm;
NagError      fail;
/*Double scalar and array declarations */
static double ruser[1] = {-1.0};
double        *resid = 0, *sigma = 0, *u = 0, *v = 0;
Nag_OrderType order;

INIT_FAIL(fail);

printf("nag_real_partial_svd (f02wgc) Example Program Results\n\n");

/* For communication with user-supplied functions: */
comm.user = ruser;

/*     Skip heading in data file*/
scanf("%*[^\n] ");
scanf("%ld%ld%ld%ld%*[^\n]", &m, &n, &nev, &ncv);

#ifdef NAG_COLUMN_MAJOR
order = Nag_ColMajor;
pdu = m;
pdv = n;
#define U(I, J) u[(J-1)*pdu + I-1]
#define V(I, J) v[(J-1)*pdv + I-1]
#else
order = Nag_RowMajor;
pdu = ncv;
pdv = ncv;
#define U(I, J) u[(I-1)*pdu + J-1]
#define V(I, J) v[(I-1)*pdv + J-1]
#endif

if (!(resid = NAG_ALLOC(m, double)) ||
!(sigma = NAG_ALLOC(ncv, double)) ||
!(u = NAG_ALLOC(m*ncv, double)) ||
!(v = NAG_ALLOC(n*ncv, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/*
* nag_real_partial_svd (f02wgc)
* Computes leading terms in the singular value decomposition of
* a real general matrix; also computes corresponding left and right
* singular vectors.
*/
nag_real_partial_svd(order, m, n, nev, ncv, av, &nconv, sigma, u, pdu,
v, pdv, resid, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_real_partial_svd (f02wgc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Print computed residuals*/
printf("%s\n", "  Singular Value    Residual");
for (i = 0; i < nconv; i++)
printf("%10.5f %16.2g\n", sigma[i], resid[i]);
printf("\n");

END:
NAG_FREE(resid);
NAG_FREE(sigma);
NAG_FREE(u);
NAG_FREE(v);

return exit_status;
}

static void NAG_CALL av(Integer *iflag, Integer m, Integer n, const double x[],
double ax[], Nag_Comm *comm)
{
Integer i, j;
double  one = 1.0, zero = 0.0;
double  h, k, s, t;

/* Matrix vector multiply subroutine Computing  w <- A*x or w <- Trans(A)*x.*/
if (comm->user[0] == -1.0)
{
printf("(User-supplied callback av, first invocation.)\n");
comm->user[0] = 0.0;
}
h = one/(double)(m+1);
k = one/(double)(n+1);
if (*iflag == 1)
{
for (i = 0; i < m; i++)
ax[i] = zero;
t = zero;
for (j = 0; j < n; j++)
{
t = t+k;
s = zero;
for (i = 0; i < MIN(j+1, m); i++)
{
s = s+h;
ax[i] = ax[i]+k*s*(t-one)*x[j];
}
for (i = j+1; i < m; i++)
{
s = s+h;
ax[i] = ax[i]+k*t*(s-one)*x[j];
}
}
}
else
{
for (i = 0; i < n; i++)
ax[i] = zero;
t = zero;
for (j = 0; j < n; j++)
{
t = t+k;
s = zero;
for (i = 0; i < MIN(j+1, m); i++)
{
s = s+h;
ax[j] = ax[j]+k*s*(t-one)*x[i];
}
for (i = j+1; i < m; i++)
{
s = s+h;
ax[j] = ax[j]+k*t*(s-one)*x[i];
}
}
}
return;
}