blender/intern/opennl/superlu/spanel_bmod.c

450 lines
13 KiB
C
Raw Normal View History

/*
* -- SuperLU routine (version 3.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* October 15, 2003
*
*/
/*
Copyright (c) 1994 by Xerox Corporation. All rights reserved.
THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
Permission is hereby granted to use or copy this program for any
purpose, provided the above notices are retained on all copies.
Permission to modify the code and to distribute modified code is
granted, provided the above notices are retained, and a notice that
the code was modified is included with the above copyright notice.
*/
#include <stdio.h>
#include <stdlib.h>
#include "ssp_defs.h"
/*
* Function prototypes
*/
void slsolve(int, int, float *, float *);
void smatvec(int, int, int, float *, float *, float *);
extern void scheck_tempv();
void
spanel_bmod (
const int m, /* in - number of rows in the matrix */
const int w, /* in */
const int jcol, /* in */
const int nseg, /* in */
float *dense, /* out, of size n by w */
float *tempv, /* working array */
int *segrep, /* in */
int *repfnz, /* in, of size n by w */
GlobalLU_t *Glu, /* modified */
SuperLUStat_t *stat /* output */
)
{
/*
* Purpose
* =======
*
* Performs numeric block updates (sup-panel) in topological order.
* It features: col-col, 2cols-col, 3cols-col, and sup-col updates.
* Special processing on the supernodal portion of L\U[*,j]
*
* Before entering this routine, the original nonzeros in the panel
* were already copied into the spa[m,w].
*
* Updated/Output parameters-
* dense[0:m-1,w]: L[*,j:j+w-1] and U[*,j:j+w-1] are returned
* collectively in the m-by-w vector dense[*].
*
*/
#ifdef USE_VENDOR_BLAS
#ifdef _CRAY
_fcd ftcs1 = _cptofcd("L", strlen("L")),
ftcs2 = _cptofcd("N", strlen("N")),
ftcs3 = _cptofcd("U", strlen("U"));
#endif
int incx = 1, incy = 1;
float alpha, beta;
#endif
register int k, ksub;
int fsupc, nsupc, nsupr, nrow;
int krep, krep_ind;
float ukj, ukj1, ukj2;
int luptr, luptr1, luptr2;
int segsze;
int block_nrow; /* no of rows in a block row */
register int lptr; /* Points to the row subscripts of a supernode */
int kfnz, irow, no_zeros;
register int isub, isub1, i;
register int jj; /* Index through each column in the panel */
int *xsup, *supno;
int *lsub, *xlsub;
float *lusup;
int *xlusup;
int *repfnz_col; /* repfnz[] for a column in the panel */
float *dense_col; /* dense[] for a column in the panel */
float *tempv1; /* Used in 1-D update */
float *TriTmp, *MatvecTmp; /* used in 2-D update */
float zero = 0.0;
register int ldaTmp;
register int r_ind, r_hi;
static int first = 1, maxsuper, rowblk, colblk;
flops_t *ops = stat->ops;
xsup = Glu->xsup;
supno = Glu->supno;
lsub = Glu->lsub;
xlsub = Glu->xlsub;
lusup = Glu->lusup;
xlusup = Glu->xlusup;
if ( first ) {
maxsuper = sp_ienv(3);
rowblk = sp_ienv(4);
colblk = sp_ienv(5);
first = 0;
}
ldaTmp = maxsuper + rowblk;
/*
* For each nonz supernode segment of U[*,j] in topological order
*/
k = nseg - 1;
for (ksub = 0; ksub < nseg; ksub++) { /* for each updating supernode */
/* krep = representative of current k-th supernode
* fsupc = first supernodal column
* nsupc = no of columns in a supernode
* nsupr = no of rows in a supernode
*/
krep = segrep[k--];
fsupc = xsup[supno[krep]];
nsupc = krep - fsupc + 1;
nsupr = xlsub[fsupc+1] - xlsub[fsupc];
nrow = nsupr - nsupc;
lptr = xlsub[fsupc];
krep_ind = lptr + nsupc - 1;
repfnz_col = repfnz;
dense_col = dense;
if ( nsupc >= colblk && nrow > rowblk ) { /* 2-D block update */
TriTmp = tempv;
/* Sequence through each column in panel -- triangular solves */
for (jj = jcol; jj < jcol + w; jj++,
repfnz_col += m, dense_col += m, TriTmp += ldaTmp ) {
kfnz = repfnz_col[krep];
if ( kfnz == EMPTY ) continue; /* Skip any zero segment */
segsze = krep - kfnz + 1;
luptr = xlusup[fsupc];
ops[TRSV] += segsze * (segsze - 1);
ops[GEMV] += 2 * nrow * segsze;
/* Case 1: Update U-segment of size 1 -- col-col update */
if ( segsze == 1 ) {
ukj = dense_col[lsub[krep_ind]];
luptr += nsupr*(nsupc-1) + nsupc;
for (i = lptr + nsupc; i < xlsub[fsupc+1]; i++) {
irow = lsub[i];
dense_col[irow] -= ukj * lusup[luptr];
++luptr;
}
} else if ( segsze <= 3 ) {
ukj = dense_col[lsub[krep_ind]];
ukj1 = dense_col[lsub[krep_ind - 1]];
luptr += nsupr*(nsupc-1) + nsupc-1;
luptr1 = luptr - nsupr;
if ( segsze == 2 ) {
ukj -= ukj1 * lusup[luptr1];
dense_col[lsub[krep_ind]] = ukj;
for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
irow = lsub[i];
luptr++; luptr1++;
dense_col[irow] -= (ukj*lusup[luptr]
+ ukj1*lusup[luptr1]);
}
} else {
ukj2 = dense_col[lsub[krep_ind - 2]];
luptr2 = luptr1 - nsupr;
ukj1 -= ukj2 * lusup[luptr2-1];
ukj = ukj - ukj1*lusup[luptr1] - ukj2*lusup[luptr2];
dense_col[lsub[krep_ind]] = ukj;
dense_col[lsub[krep_ind-1]] = ukj1;
for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
irow = lsub[i];
luptr++; luptr1++; luptr2++;
dense_col[irow] -= ( ukj*lusup[luptr]
+ ukj1*lusup[luptr1] + ukj2*lusup[luptr2] );
}
}
} else { /* segsze >= 4 */
/* Copy U[*,j] segment from dense[*] to TriTmp[*], which
holds the result of triangular solves. */
no_zeros = kfnz - fsupc;
isub = lptr + no_zeros;
for (i = 0; i < segsze; ++i) {
irow = lsub[isub];
TriTmp[i] = dense_col[irow]; /* Gather */
++isub;
}
/* start effective triangle */
luptr += nsupr * no_zeros + no_zeros;
#ifdef USE_VENDOR_BLAS
#ifdef _CRAY
STRSV( ftcs1, ftcs2, ftcs3, &segsze, &lusup[luptr],
&nsupr, TriTmp, &incx );
#else
strsv_( "L", "N", "U", &segsze, &lusup[luptr],
&nsupr, TriTmp, &incx );
#endif
#else
slsolve ( nsupr, segsze, &lusup[luptr], TriTmp );
#endif
} /* else ... */
} /* for jj ... end tri-solves */
/* Block row updates; push all the way into dense[*] block */
for ( r_ind = 0; r_ind < nrow; r_ind += rowblk ) {
r_hi = SUPERLU_MIN(nrow, r_ind + rowblk);
block_nrow = SUPERLU_MIN(rowblk, r_hi - r_ind);
luptr = xlusup[fsupc] + nsupc + r_ind;
isub1 = lptr + nsupc + r_ind;
repfnz_col = repfnz;
TriTmp = tempv;
dense_col = dense;
/* Sequence through each column in panel -- matrix-vector */
for (jj = jcol; jj < jcol + w; jj++,
repfnz_col += m, dense_col += m, TriTmp += ldaTmp) {
kfnz = repfnz_col[krep];
if ( kfnz == EMPTY ) continue; /* Skip any zero segment */
segsze = krep - kfnz + 1;
if ( segsze <= 3 ) continue; /* skip unrolled cases */
/* Perform a block update, and scatter the result of
matrix-vector to dense[]. */
no_zeros = kfnz - fsupc;
luptr1 = luptr + nsupr * no_zeros;
MatvecTmp = &TriTmp[maxsuper];
#ifdef USE_VENDOR_BLAS
alpha = one;
beta = zero;
#ifdef _CRAY
SGEMV(ftcs2, &block_nrow, &segsze, &alpha, &lusup[luptr1],
&nsupr, TriTmp, &incx, &beta, MatvecTmp, &incy);
#else
sgemv_("N", &block_nrow, &segsze, &alpha, &lusup[luptr1],
&nsupr, TriTmp, &incx, &beta, MatvecTmp, &incy);
#endif
#else
smatvec(nsupr, block_nrow, segsze, &lusup[luptr1],
TriTmp, MatvecTmp);
#endif
/* Scatter MatvecTmp[*] into SPA dense[*] temporarily
* such that MatvecTmp[*] can be re-used for the
* the next blok row update. dense[] will be copied into
* global store after the whole panel has been finished.
*/
isub = isub1;
for (i = 0; i < block_nrow; i++) {
irow = lsub[isub];
dense_col[irow] -= MatvecTmp[i];
MatvecTmp[i] = zero;
++isub;
}
} /* for jj ... */
} /* for each block row ... */
/* Scatter the triangular solves into SPA dense[*] */
repfnz_col = repfnz;
TriTmp = tempv;
dense_col = dense;
for (jj = jcol; jj < jcol + w; jj++,
repfnz_col += m, dense_col += m, TriTmp += ldaTmp) {
kfnz = repfnz_col[krep];
if ( kfnz == EMPTY ) continue; /* Skip any zero segment */
segsze = krep - kfnz + 1;
if ( segsze <= 3 ) continue; /* skip unrolled cases */
no_zeros = kfnz - fsupc;
isub = lptr + no_zeros;
for (i = 0; i < segsze; i++) {
irow = lsub[isub];
dense_col[irow] = TriTmp[i];
TriTmp[i] = zero;
++isub;
}
} /* for jj ... */
} else { /* 1-D block modification */
/* Sequence through each column in the panel */
for (jj = jcol; jj < jcol + w; jj++,
repfnz_col += m, dense_col += m) {
kfnz = repfnz_col[krep];
if ( kfnz == EMPTY ) continue; /* Skip any zero segment */
segsze = krep - kfnz + 1;
luptr = xlusup[fsupc];
ops[TRSV] += segsze * (segsze - 1);
ops[GEMV] += 2 * nrow * segsze;
/* Case 1: Update U-segment of size 1 -- col-col update */
if ( segsze == 1 ) {
ukj = dense_col[lsub[krep_ind]];
luptr += nsupr*(nsupc-1) + nsupc;
for (i = lptr + nsupc; i < xlsub[fsupc+1]; i++) {
irow = lsub[i];
dense_col[irow] -= ukj * lusup[luptr];
++luptr;
}
} else if ( segsze <= 3 ) {
ukj = dense_col[lsub[krep_ind]];
luptr += nsupr*(nsupc-1) + nsupc-1;
ukj1 = dense_col[lsub[krep_ind - 1]];
luptr1 = luptr - nsupr;
if ( segsze == 2 ) {
ukj -= ukj1 * lusup[luptr1];
dense_col[lsub[krep_ind]] = ukj;
for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
irow = lsub[i];
++luptr; ++luptr1;
dense_col[irow] -= (ukj*lusup[luptr]
+ ukj1*lusup[luptr1]);
}
} else {
ukj2 = dense_col[lsub[krep_ind - 2]];
luptr2 = luptr1 - nsupr;
ukj1 -= ukj2 * lusup[luptr2-1];
ukj = ukj - ukj1*lusup[luptr1] - ukj2*lusup[luptr2];
dense_col[lsub[krep_ind]] = ukj;
dense_col[lsub[krep_ind-1]] = ukj1;
for (i = lptr + nsupc; i < xlsub[fsupc+1]; ++i) {
irow = lsub[i];
++luptr; ++luptr1; ++luptr2;
dense_col[irow] -= ( ukj*lusup[luptr]
+ ukj1*lusup[luptr1] + ukj2*lusup[luptr2] );
}
}
} else { /* segsze >= 4 */
/*
* Perform a triangular solve and block update,
* then scatter the result of sup-col update to dense[].
*/
no_zeros = kfnz - fsupc;
/* Copy U[*,j] segment from dense[*] to tempv[*]:
* The result of triangular solve is in tempv[*];
* The result of matrix vector update is in dense_col[*]
*/
isub = lptr + no_zeros;
for (i = 0; i < segsze; ++i) {
irow = lsub[isub];
tempv[i] = dense_col[irow]; /* Gather */
++isub;
}
/* start effective triangle */
luptr += nsupr * no_zeros + no_zeros;
#ifdef USE_VENDOR_BLAS
#ifdef _CRAY
STRSV( ftcs1, ftcs2, ftcs3, &segsze, &lusup[luptr],
&nsupr, tempv, &incx );
#else
strsv_( "L", "N", "U", &segsze, &lusup[luptr],
&nsupr, tempv, &incx );
#endif
luptr += segsze; /* Dense matrix-vector */
tempv1 = &tempv[segsze];
alpha = one;
beta = zero;
#ifdef _CRAY
SGEMV( ftcs2, &nrow, &segsze, &alpha, &lusup[luptr],
&nsupr, tempv, &incx, &beta, tempv1, &incy );
#else
sgemv_( "N", &nrow, &segsze, &alpha, &lusup[luptr],
&nsupr, tempv, &incx, &beta, tempv1, &incy );
#endif
#else
slsolve ( nsupr, segsze, &lusup[luptr], tempv );
luptr += segsze; /* Dense matrix-vector */
tempv1 = &tempv[segsze];
smatvec (nsupr, nrow, segsze, &lusup[luptr], tempv, tempv1);
#endif
/* Scatter tempv[*] into SPA dense[*] temporarily, such
* that tempv[*] can be used for the triangular solve of
* the next column of the panel. They will be copied into
* ucol[*] after the whole panel has been finished.
*/
isub = lptr + no_zeros;
for (i = 0; i < segsze; i++) {
irow = lsub[isub];
dense_col[irow] = tempv[i];
tempv[i] = zero;
isub++;
}
/* Scatter the update from tempv1[*] into SPA dense[*] */
/* Start dense rectangular L */
for (i = 0; i < nrow; i++) {
irow = lsub[isub];
dense_col[irow] -= tempv1[i];
tempv1[i] = zero;
++isub;
}
} /* else segsze>=4 ... */
} /* for each column in the panel... */
} /* else 1-D update ... */
} /* for each updating supernode ... */
}