本文整理汇总了C++中saxpy函数的典型用法代码示例。如果您正苦于以下问题:C++ saxpy函数的具体用法?C++ saxpy怎么用?C++ saxpy使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了saxpy函数的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: axpy
static vl::Error
axpy(vl::Context & context,
ptrdiff_t n,
type alpha,
type const *x, ptrdiff_t incx,
type *y, ptrdiff_t incy)
{
saxpy(&n,
&alpha,
(float*)x, &incx,
(float*)y, &incy) ;
return vl::vlSuccess ;
}
示例2: MocCompute2WayPartitionParams
/*************************************************************************
* This function computes the initial id/ed
**************************************************************************/
void MocCompute2WayPartitionParams(CtrlType *ctrl, GraphType *graph)
{
int i, j, /*k, l,*/ nvtxs, ncon, nbnd, mincut;
idxtype *xadj, *adjncy, *adjwgt;
float *nvwgt, *npwgts;
idxtype *id, *ed, *where;
idxtype *bndptr, *bndind;
int me/*, other*/;
nvtxs = graph->nvtxs;
ncon = graph->ncon;
xadj = graph->xadj;
nvwgt = graph->nvwgt;
adjncy = graph->adjncy;
adjwgt = graph->adjwgt;
where = graph->where;
npwgts = sset(2*ncon, 0.0, graph->npwgts);
id = idxset(nvtxs, 0, graph->id);
ed = idxset(nvtxs, 0, graph->ed);
bndptr = idxset(nvtxs, -1, graph->bndptr);
bndind = graph->bndind;
/*------------------------------------------------------------
/ Compute now the id/ed degrees
/------------------------------------------------------------*/
nbnd = mincut = 0;
for (i=0; i<nvtxs; i++) {
ASSERT(where[i] >= 0 && where[i] <= 1);
me = where[i];
saxpy(ncon, 1.0, nvwgt+i*ncon, 1, npwgts+me*ncon, 1);
for (j=xadj[i]; j<xadj[i+1]; j++) {
if (me == where[adjncy[j]])
id[i] += adjwgt[j];
else
ed[i] += adjwgt[j];
}
if (ed[i] > 0 || xadj[i] == xadj[i+1]) {
mincut += ed[i];
bndptr[i] = nbnd;
bndind[nbnd++] = i;
}
}
graph->mincut = mincut/2;
graph->nbnd = nbnd;
}
示例3: interpolate_initial
//! @{
virtual void interpolate_initial(shared_ptr<ISweeper<time>> dst,
shared_ptr<const ISweeper<time>> src) override
{
auto& fine = as_encap_sweeper(dst);
auto& crse = as_encap_sweeper(src);
auto crse_factory = crse.get_factory();
auto fine_factory = fine.get_factory();
auto crse_delta = crse_factory->create(solution);
this->restrict(crse_delta, fine.get_start_state());
crse_delta->saxpy(-1.0, crse.get_start_state());
auto fine_delta = fine_factory->create(solution);
this->interpolate(fine_delta, crse_delta);
fine.get_start_state()->saxpy(-1.0, fine_delta);
fine.reevaluate(true);
}
示例4: main
int main()
{
uint32_t num_cpus = std::thread::hardware_concurrency();
std::vector<std::thread>threads(num_cpus);
float *S = new float[num_cpus], a = 1.0f, b = 2.0f, *X = new float[num_cpus], *Y = new float[num_cpus];
for(uint32_t i=0;i<num_cpus;i++){
X[i] = 1.0f*i;
Y[i] = 2.0f*(i+1);
Saxpy saxpy(std::ref(S[i]), a, X[i], b, Y[i]);
threads[i] = std::thread(saxpy);
}
for(uint32_t i=0;i<num_cpus;i++){
threads[i].join();
std::cout<<S[i]<<std::endl;
}
}
示例5: main
int main (int argc, char * argv[])
{
const int N = 16;
const int iters = 1;
float *x, *y, *z;
posix_memalign((void **)&x, VECTOR_SIZE, N*sizeof(float));
posix_memalign((void **)&y, VECTOR_SIZE, N*sizeof(float));
posix_memalign((void **)&z, VECTOR_SIZE, N*sizeof(float));
float a = 0.93f;
int i, j;
for (i=0; i<N; i++)
{
x[i] = i+1;
y[i] = i-1;
z[i] = 0.0f;
}
for (i=0; i<iters; i++)
{
saxpy(x, y, z, a, N);
}
for (i=0; i<N; i++)
{
if (z[i] != (a * x[i] + y[i]))
{
printf("Error\n");
return (1);
}
}
printf("SUCCESS!\n");
return 0;
}
示例6: main
int main(int argc, char* argv[]) {
float a, x[N], y[N];
a=5;
int i;
for (i=0; i<N; ++i){
x[i]=i;
y[i]=i+2;
}
saxpy(N, a, x, y);
#pragma omp taskwait
//Check results
for (i=0; i<N; ++i){
if (y[i]!=a*i+(i+2)){
printf("Error when checking results, in position %d\n",i);
return -1;
}
}
printf("Results are correct\n");
return 0;
}
示例7: interpolate
virtual void interpolate(shared_ptr<ISweeper<time>> dst,
shared_ptr<const ISweeper<time>> src,
bool interp_initial) override
{
auto& fine = as_encap_sweeper(dst);
auto& crse = as_encap_sweeper(src);
if (tmat.rows() == 0) {
tmat = pfasst::quadrature::compute_interp<time>(fine.get_nodes(), crse.get_nodes());
}
if (interp_initial) {
this->interpolate_initial(dst, src);
}
size_t nfine = fine.get_nodes().size();
size_t ncrse = crse.get_nodes().size();
auto crse_factory = crse.get_factory();
auto fine_factory = fine.get_factory();
EncapVecT fine_state(nfine), fine_delta(ncrse);
for (size_t m = 0; m < nfine; m++) { fine_state[m] = fine.get_state(m); }
for (size_t m = 0; m < ncrse; m++) { fine_delta[m] = fine_factory->create(solution); }
auto crse_delta = crse_factory->create(solution);
for (size_t m = 0; m < ncrse; m++) {
crse_delta->copy(crse.get_state(m));
crse_delta->saxpy(-1.0, crse.get_saved_state(m));
interpolate(fine_delta[m], crse_delta);
}
fine.get_state(0)->mat_apply(fine_state, 1.0, tmat, fine_delta, false);
fine.reevaluate();
}
示例8: operator
void operator()(float a, float *x, float *y, std::size_t n)
{
saxpy(a,x,y,n);
}
示例9: sgefa
int sgefa ( float a[], int lda, int n, int ipvt[] )
/*******************************************************************************/
/*
Purpose:
SGEFA factors a matrix by gaussian elimination.
Discussion:
Matrix references which would, mathematically, be written A(I,J)
must be written here as:
* A[I+J*LDA], when the value is needed, or
* A+I+J*LDA, when the address is needed.
Modified:
07 March 2008
Author:
FORTRAN77 original version by Cleve Moler.
C version by John Burkardt.
Reference:
Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
LINPACK User's Guide,
SIAM, 1979,
ISBN13: 978-0-898711-72-1,
LC: QA214.L56.
Parameters:
Input/output, float A[LDA*N]. On input, the matrix to be factored.
On output, an upper triangular matrix and the multipliers which were
used to obtain it. The factorization can be written A = L * U where
L is a product of permutation and unit lower triangular matrices and
U is upper triangular.
Input, int LDA, the leading dimension of the matrix.
Input, int N, the order of the matrix.
Output, int IPVT[N], the pivot indices.
Output, int SGEFA, indicates singularity.
If 0, this is the normal value, and the algorithm succeeded.
If K, then on the K-th elimination step, a zero pivot was encountered.
The matrix is numerically not invertible.
*/
{
int j;
int info;
int k;
int kp1;
int l;
int nm1;
float t;
info = 0;
for ( k = 1; k <= n - 1; k++ )
{
/*
Find l = pivot index.
*/
l = isamax ( n-k+1, &a[k-1+(k-1)*lda], 1 ) + k - 1;
ipvt[k-1] = l;
/*
Zero pivot implies this column already triangularized.
*/
if ( a[l-1+(k-1)*lda] != 0.0 )
{
/*
Interchange if necessary.
*/
if ( l != k )
{
t = a[l-1+(k-1)*lda];
a[l-1+(k-1)*lda] = a[k-1+(k-1)*lda];
a[k-1+(k-1)*lda] = t;
}
/*
Compute multipliers.
*/
t = - 1.0 / a[k-1+(k-1)*lda];
sscal ( n-k, t, &a[k+(k-1)*lda], 1 );
/*
Row elimination with column indexing.
*/
for ( j = k + 1; j <= n; j++ )
{
t = a[l-1+(j-1)*lda];
if (l != k)
{
a[l-1+(j-1)*lda] = a[k-1+(j-1)*lda];
a[k-1+(j-1)*lda] = t;
}
saxpy ( n-k, t, &a[k+(k-1)*lda], 1, &a[k+(j-1)*lda], 1 );
//.........这里部分代码省略.........
示例10: main
int
main()
{
#define N 8
int i;
float x_ref[N], y_ref[N];
float x[N], y[N];
cublasHandle_t h;
float a = 2.0;
for (i = 0; i < N; i++)
{
x[i] = x_ref[i] = 4.0 + i;
y[i] = y_ref[i] = 3.0;
}
saxpy (N, a, x_ref, y_ref);
cublasCreate (&h);
#pragma acc data copyin (x[0:N]) copy (y[0:N])
{
#pragma acc host_data use_device (x, y)
{
cublasSaxpy (h, N, &a, x, 1, y, 1);
}
}
validate_results (N, y, y_ref);
#pragma acc data create (x[0:N]) copyout (y[0:N])
{
#pragma acc kernels
for (i = 0; i < N; i++)
y[i] = 3.0;
#pragma acc host_data use_device (x, y)
{
cublasSaxpy (h, N, &a, x, 1, y, 1);
}
}
cublasDestroy (h);
validate_results (N, y, y_ref);
for (i = 0; i < N; i++)
y[i] = 3.0;
/* There's no need to use host_data here. */
#pragma acc data copyin (x[0:N]) copyin (a) copy (y[0:N])
{
#pragma acc parallel present (x[0:N]) pcopy (y[0:N]) present (a)
saxpy (N, a, x, y);
}
validate_results (N, y, y_ref);
/* Exercise host_data with data transferred with acc enter data. */
for (i = 0; i < N; i++)
y[i] = 3.0;
#pragma acc enter data copyin (x, a, y)
#pragma acc parallel present (x[0:N]) pcopy (y[0:N]) present (a)
{
saxpy (N, a, x, y);
}
#pragma acc exit data delete (x, a) copyout (y)
validate_results (N, y, y_ref);
return 0;
}
示例11: CreateCoarseGraphNoMask
/*************************************************************************
* This function creates the coarser graph
**************************************************************************/
void CreateCoarseGraphNoMask(CtrlType *ctrl, GraphType *graph, int cnvtxs, idxtype *match, idxtype *perm)
{
int i, j, k, m, istart, iend, nvtxs, nedges, ncon, cnedges, v, u, dovsize;
idxtype *xadj, *vwgt, *vsize, *adjncy, *adjwgt, *adjwgtsum, *auxadj;
idxtype *cmap, *htable;
idxtype *cxadj, *cvwgt, *cvsize, *cadjncy, *cadjwgt, *cadjwgtsum;
float *nvwgt, *cnvwgt;
GraphType *cgraph;
dovsize = (ctrl->optype == OP_KVMETIS ? 1 : 0);
IFSET(ctrl->dbglvl, DBG_TIME, starttimer(ctrl->ContractTmr));
nvtxs = graph->nvtxs;
ncon = graph->ncon;
xadj = graph->xadj;
vwgt = graph->vwgt;
vsize = graph->vsize;
nvwgt = graph->nvwgt;
adjncy = graph->adjncy;
adjwgt = graph->adjwgt;
adjwgtsum = graph->adjwgtsum;
cmap = graph->cmap;
/* Initialize the coarser graph */
cgraph = SetUpCoarseGraph(graph, cnvtxs, dovsize);
cxadj = cgraph->xadj;
cvwgt = cgraph->vwgt;
cvsize = cgraph->vsize;
cnvwgt = cgraph->nvwgt;
cadjwgtsum = cgraph->adjwgtsum;
cadjncy = cgraph->adjncy;
cadjwgt = cgraph->adjwgt;
htable = idxset(cnvtxs, -1, idxwspacemalloc(ctrl, cnvtxs));
iend = xadj[nvtxs];
auxadj = ctrl->wspace.auxcore;
memcpy(auxadj, adjncy, iend*sizeof(idxtype));
for (i=0; i<iend; i++)
auxadj[i] = cmap[auxadj[i]];
cxadj[0] = cnvtxs = cnedges = 0;
for (i=0; i<nvtxs; i++) {
v = perm[i];
if (cmap[v] != cnvtxs)
continue;
u = match[v];
if (ncon == 1)
cvwgt[cnvtxs] = vwgt[v];
else
scopy(ncon, nvwgt+v*ncon, cnvwgt+cnvtxs*ncon);
if (dovsize)
cvsize[cnvtxs] = vsize[v];
cadjwgtsum[cnvtxs] = adjwgtsum[v];
nedges = 0;
istart = xadj[v];
iend = xadj[v+1];
for (j=istart; j<iend; j++) {
k = auxadj[j];
if ((m = htable[k]) == -1) {
cadjncy[nedges] = k;
cadjwgt[nedges] = adjwgt[j];
htable[k] = nedges++;
}
else {
cadjwgt[m] += adjwgt[j];
}
}
if (v != u) {
if (ncon == 1)
cvwgt[cnvtxs] += vwgt[u];
else
saxpy(ncon, 1.0, nvwgt+u*ncon, 1, cnvwgt+cnvtxs*ncon, 1);
if (dovsize)
cvsize[cnvtxs] += vsize[u];
cadjwgtsum[cnvtxs] += adjwgtsum[u];
istart = xadj[u];
iend = xadj[u+1];
for (j=istart; j<iend; j++) {
k = auxadj[j];
if ((m = htable[k]) == -1) {
cadjncy[nedges] = k;
cadjwgt[nedges] = adjwgt[j];
htable[k] = nedges++;
}
else {
//.........这里部分代码省略.........
示例12: MCRandom_KWayEdgeRefineHorizontal
//.........这里部分代码省略.........
if (gain >= 0 &&
(AreAllHVwgtsBelow(ncon, 1.0, npwgts+to*ncon, 1.0, nvwgt, maxwgt+to*ncon) ||
IsHBalanceBetterFT(ncon, nparts, npwgts+from*ncon, npwgts+to*ncon, nvwgt, ubvec)))
break;
}
if (k == myndegrees)
continue; /* break out if you did not find a candidate */
for (j=k+1; j<myndegrees; j++) {
to = myedegrees[j].pid;
if ((myedegrees[j].ed > myedegrees[k].ed &&
(AreAllHVwgtsBelow(ncon, 1.0, npwgts+to*ncon, 1.0, nvwgt, maxwgt+to*ncon) ||
IsHBalanceBetterFT(ncon, nparts, npwgts+from*ncon, npwgts+to*ncon, nvwgt, ubvec))) ||
(myedegrees[j].ed == myedegrees[k].ed &&
IsHBalanceBetterTT(ncon, nparts, npwgts+myedegrees[k].pid*ncon, npwgts+to*ncon, nvwgt, ubvec)))
k = j;
}
to = myedegrees[k].pid;
if (myedegrees[k].ed-myrinfo->id == 0
&& !IsHBalanceBetterFT(ncon, nparts, npwgts+from*ncon, npwgts+to*ncon, nvwgt, ubvec)
&& AreAllHVwgtsBelow(ncon, 1.0, npwgts+from*ncon, 0.0, npwgts+from*ncon, maxwgt+from*ncon))
continue;
/*=====================================================================
* If we got here, we can now move the vertex from 'from' to 'to'
*======================================================================*/
graph->mincut -= myedegrees[k].ed-myrinfo->id;
IFSET(ctrl->dbglvl, DBG_MOVEINFO, printf("\t\tMoving %6d to %3d. Gain: %4d. Cut: %6d\n", i, to, myedegrees[k].ed-myrinfo->id, graph->mincut));
/* Update where, weight, and ID/ED information of the vertex you moved */
saxpy(ncon, 1.0, nvwgt, 1, npwgts+to*ncon, 1);
saxpy(ncon, -1.0, nvwgt, 1, npwgts+from*ncon, 1);
where[i] = to;
myrinfo->ed += myrinfo->id-myedegrees[k].ed;
SWAP(myrinfo->id, myedegrees[k].ed, j);
if (myedegrees[k].ed == 0)
myedegrees[k] = myedegrees[--myrinfo->ndegrees];
else
myedegrees[k].pid = from;
if (myrinfo->ed-myrinfo->id < 0)
BNDDelete(nbnd, bndind, bndptr, i);
/* Update the degrees of adjacent vertices */
for (j=xadj[i]; j<xadj[i+1]; j++) {
ii = adjncy[j];
me = where[ii];
myrinfo = graph->rinfo+ii;
if (myrinfo->edegrees == NULL) {
myrinfo->edegrees = ctrl->wspace.edegrees+ctrl->wspace.cdegree;
ctrl->wspace.cdegree += xadj[ii+1]-xadj[ii];
}
myedegrees = myrinfo->edegrees;
ASSERT(CheckRInfo(myrinfo));
if (me == from) {
INC_DEC(myrinfo->ed, myrinfo->id, adjwgt[j]);
if (myrinfo->ed-myrinfo->id >= 0 && bndptr[ii] == -1)
BNDInsert(nbnd, bndind, bndptr, ii);
}
示例13: sgesl
void sgesl ( float a[], int lda, int n, int ipvt[], float b[], int job )
/******************************************************************************/
/*
Purpose:
SGESL solves a real general linear system A * X = B.
Discussion:
SGESL can solve either of the systems A * X = B or A' * X = B.
The system matrix must have been factored by SGECO or SGEFA.
A division by zero will occur if the input factor contains a
zero on the diagonal. Technically this indicates singularity
but it is often caused by improper arguments or improper
setting of LDA. It will not occur if the subroutines are
called correctly and if SGECO has set 0.0 < RCOND
or SGEFA has set INFO == 0.
Modified:
04 April 2006
Author:
FORTRAN77 original by Dongarra, Moler, Bunch and Stewart.
C translation by John Burkardt.
Reference:
Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
LINPACK User's Guide,
SIAM, (Society for Industrial and Applied Mathematics),
3600 University City Science Center,
Philadelphia, PA, 19104-2688.
ISBN: 0-89871-172-X
Parameters:
Input, float A[LDA*N], the output from SGECO or SGEFA.
Input, int LDA, the leading dimension of A.
Input, int N, the order of the matrix A.
Input, int IPVT[N], the pivot vector from SGECO or SGEFA.
Input/output, float B[N].
On input, the right hand side vector.
On output, the solution vector.
Input, int JOB.
0, solve A * X = B;
nonzero, solve A' * X = B.
*/
{
int k;
int l;
float t;
/*
Solve A * X = B.
*/
if ( job == 0 )
{
for ( k = 1; k <= n-1; k++ )
{
l = ipvt[k-1];
t = b[l-1];
if ( l != k )
{
b[l-1] = b[k-1];
b[k-1] = t;
}
saxpy ( n-k, t, a+k+(k-1)*lda, 1, b+k, 1 );
}
for ( k = n; 1 <= k; k-- )
{
b[k-1] = b[k-1] / a[k-1+(k-1)*lda];
t = -b[k-1];
saxpy ( k-1, t, a+0+(k-1)*lda, 1, b, 1 );
}
}
/*
Solve A' * X = B.
*/
else
{
for ( k = 1; k <= n; k++ )
{
t = sdot ( k-1, a+0+(k-1)*lda, 1, b, 1 );
b[k-1] = ( b[k-1] - t ) / a[k-1+(k-1)*lda];
}
for ( k = n-1; 1 <= k; k-- )
{
b[k-1] = b[k-1] + sdot ( n-k, a+k+(k-1)*lda, 1, b+k, 1 );
//.........这里部分代码省略.........
示例14: ConjGrad2
/*************************************************************************
* This function implements the CG solver used during the directed diffusion
**************************************************************************/
void ConjGrad2(MatrixType *A, floattype *b, floattype *x, floattype tol, floattype *workspace)
{
int i, k, n;
floattype *p, *r, *q, *z, *M;
floattype alpha, beta, rho, rho_1 = -1.0, error, bnrm2, tmp;
idxtype *rowptr, *colind;
floattype *values;
n = A->nrows;
rowptr = A->rowptr;
colind = A->colind;
values = A->values;
/* Initial Setup */
p = workspace;
r = workspace + n;
q = workspace + 2*n;
z = workspace + 3*n;
M = workspace + 4*n;
for (i=0; i<n; i++) {
x[i] = 0.0;
if (values[rowptr[i]] != 0.0)
M[i] = 1.0/values[rowptr[i]];
else
M[i] = 0.0;
}
/* r = b - Ax */
mvMult2(A, x, r);
for (i=0; i<n; i++)
r[i] = b[i]-r[i];
bnrm2 = snorm2(n, b);
if (bnrm2 > 0.0) {
error = snorm2(n, r) / bnrm2;
if (error > tol) {
/* Begin Iterations */
for (k=0; k<n; k++) {
for (i=0; i<n; i++)
z[i] = r[i]*M[i];
rho = sdot(n, r, z);
if (k == 0)
scopy(n, z, p);
else {
if (rho_1 != 0.0)
beta = rho/rho_1;
else
beta = 0.0;
for (i=0; i<n; i++)
p[i] = z[i] + beta*p[i];
}
mvMult2(A, p, q); /* q = A*p */
tmp = sdot(n, p, q);
if (tmp != 0.0)
alpha = rho/tmp;
else
alpha = 0.0;
saxpy(n, alpha, p, x); /* x = x + alpha*p */
saxpy(n, -alpha, q, r); /* r = r - alpha*q */
error = snorm2(n, r) / bnrm2;
if (error < tol)
break;
rho_1 = rho;
}
}
}
}
示例15: sqrdc
//.........这里部分代码省略.........
jpvt[j] = jpvt[pl];
jpvt[pl] = j;
pl++;
}
}
pu = p-1;
for (j=p-1; j>=0; j--) {
if (jpvt[j]<0) {
jpvt[j] = -jpvt[j];
if (j!=pu) {
sswap(n,x[pu],1,x[j],1);
jp = jpvt[pu];
jpvt[pu] = jpvt[j];
jpvt[j] = jp;
}
pu--;
}
}
}
/* compute the norms of the free columns */
for (j=pl; j<=pu; j++) {
qraux[j] = snrm2(n,x[j],1);
work[j] = qraux[j];
}
/* perform the Householder reduction of x */
lup = MIN(n,p);
for (l=0; l<lup; l++) {
if (l>=pl && l<pu) {
/*
* locate the column of largest norm and
* bring it into pivot position.
*/
maxnrm = 0.0;
maxj = l;
for (j=l; j<=pu; j++) {
if (qraux[j]>maxnrm) {
maxnrm = qraux[j];
maxj = j;
}
}
if (maxj!=l) {
sswap(n,x[l],1,x[maxj],1);
qraux[maxj] = qraux[l];
work[maxj] = work[l];
jp = jpvt[maxj];
jpvt[maxj] = jpvt[l];
jpvt[l] = jp;
}
}
qraux[l] = 0.0;
if (l!=n-1) {
/*
* compute the Householder transformation
* for column l
*/
nrmxl = snrm2(n-l,&x[l][l],1);
if (nrmxl!=0.0) {
if (x[l][l]!=0.0)
nrmxl = (x[l][l]>0.0) ?
ABS(nrmxl) :
-ABS(nrmxl);
sscal(n-l,1.0/nrmxl,&x[l][l],1);
x[l][l] += 1.0;
/*
* apply the transformation to the remaining
* columns, updating the norms
*/
for (j=l+1; j<p; j++) {
t = -sdot(n-l,&x[l][l],1,&x[j][l],1)/
x[l][l];
saxpy(n-l,t,&x[l][l],1,&x[j][l],1);
if (j>=pl && j<=pu && qraux[j]!=0.0) {
tt = ABS(x[j][l])/qraux[j];
tt = 1.0-tt*tt;
tt = MAX(tt,0.0);
t = tt;
ttt = qraux[j]/work[j];
tt = 1.0+0.05*tt*ttt*ttt;
if (tt!=1.0) {
qraux[j] *= sqrt(t);
} else {
qraux[j] = snrm2(n-l-1,
&x[j][l+1],1);
work[j] = qraux[j];
}
}
}
/* save the transformation */
qraux[l] = x[l][l];
x[l][l] = -nrmxl;
}
}
}
}