From 5bb72c97b02e1cea0dbf008ec16029491956f971 Mon Sep 17 00:00:00 2001 From: jofret Date: Tue, 28 Apr 2009 06:57:25 +0000 Subject: Moving source code --- src/signalProcessing/fft/dfftmx.c | 1211 ------------------------------------- 1 file changed, 1211 deletions(-) delete mode 100644 src/signalProcessing/fft/dfftmx.c (limited to 'src/signalProcessing/fft/dfftmx.c') diff --git a/src/signalProcessing/fft/dfftmx.c b/src/signalProcessing/fft/dfftmx.c deleted file mode 100644 index f7d6ce30..00000000 --- a/src/signalProcessing/fft/dfftmx.c +++ /dev/null @@ -1,1211 +0,0 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2008-2008 - INRIA - Allan SIMON - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -#include -#include -#include "max.h" -#include "min.h" -#include "fft_internal.h" - - -/*c'est moche je sais */ - -static double* a ; -static double* b ; -static int ntot ; -static int n ; -static int nspan ; -static int isn ; -static int m ; -static int kt ; -static double* wt ; -static double* ck ; -static double* bt ; -static double* sk ; -static int* np ; -static int* nfac ; -static int inc ; -static int nt ; -static int ks ; - -static double rad ; -static double c72 ; -static double s72 ; -static double s120 ; - -static double aa ; -static double ak ; -static double akm ; -static double akp ; -static double aj ; -static double ajp ; -static double ajm ; - -static double bb ; -static double bk ; -static double bkm ; -static double bkp ; -static double bj ; -static double bjp ; -static double bjm ; -static double dr ; - -static double cd ; -static double c1 ; -static double c2 ; -static double c3 ; - -static double sd ; -static double s1 ; -static double s2 ; -static double s3 ; - -static int kspan ; -static int nn ; -static int jc ; -static int klim ; -static int jf ; -static int maxf ; -static int kk ; -static int k ; -static int k1 ; -static int k2 ; -static int k3 ; -static int k4 ; -static int mm ; -static int kspnn ; -static int i ; -static int j ; -static int jj; - -/* Prototypes */ - -static void preliminaryWork (void); -static void permute_stage1 (void); -static void permute_stage2 (void); -static void f4t_150 (void); -static void factorOf3Transform (void) ; -static void factorOf5Transform (void) ; -static void preFOtherTransform (void); -static void factorOfOtherTransform (void); -static void pre_sqFactor2NormlOrder (void); -static void nonSqFactor2NormOrder (void) ; -static void detPermutCycles (void); -static void reorderMatrix (void ) ; - -static int f4t_170 (void); -static int factorTransform (void); -static int pre_fOf2Trans (void); -static int factorOf2Transform (void); -static int factorOf4Transform (void); -static int mulByRotationFactor (void ); -static int post_sqFactor2NormlOrder (void); -static void single_sqFactor2NormlOrder (void); -static int multi_sqFactor2NormlOrder (void); - -/* End Prototypes */ - -/*note on this code all numbers alone in comment is - a reference to the corresponding goto in the original fotran code */ - -int dfftmx ( double* _pdblA , double* _pdblB , int _iNtot, int _iN, int _iNspan, - int _iIsn, int _iM, int _iKt, double* _pdblWt, double* _pdblCk, - double* _pdblBt, double* _pdblSk, int* _piNp, int* _piNfac) -{ - int retVal = 0 ; - a = _pdblA ; - b = _pdblB ; - - ntot = _iNtot ; - n = _iN ; - nspan= _iNspan ; - isn = _iIsn; - m = _iM ; - kt = _iKt ; - wt = _pdblWt ; - ck = _pdblCk; - bt = _pdblBt; - sk = _pdblSk; - np = _piNp; - nfac = _piNfac; - - inc = abs ( isn ) ; - nt = inc*ntot ; - ks = inc*nspan; - rad = atan((double)1); - c72 = cos (rad/0.6250); - s72 = sin (rad/0.6250); - s120= sqrt(0.750); - - - - preliminaryWork() ; - - - while ( retVal == 0 ) retVal = factorTransform ( ) ; - - np[0] = ks ; - - if ( kt != 0) permute_stage1 ( ) ; - - - if ( 2*kt + 1 < m ) permute_stage2 ( ) ; - - _pdblA = a ; - _pdblB = b ; - - return 0 ; -} - -/** ************************************** -Sous-Fonctions -******************************************/ - - - -/* this function only set the value of variable */ -static void preliminaryWork (void) -{ - s72 = -s72 ; - s120= -s120; - rad = -rad ; - - kspan = ks ; - nn = nt -inc ; - jc = ks/n ; - - /* sin , cos values are re-initialized each lim steps */ - - klim = 32*jc ; - i = 0; - jf = 0 ; - - if ( kt > 0 ) maxf = max ( nfac[kt-1] , nfac[m-kt-1] ); - else maxf = nfac[m-kt-1] ; - -} - - - -/*40*/ -/* this function is call as many time as dfftbi has determined factor for the size of the input vector - each time we call a transform function for each kind of factor , we begin by the smallest - factor are stored in nfac - */ - -static int factorTransform (void) -{ - - int retVal = 42; - - dr = 8 * (double)jc/(double)kspan ; - cd = 2 * sin(0.5*dr*rad)*sin(0.5*dr*rad); - sd = sin(dr*rad) ; - kk = 1 ; - i++ ; - - - - -switch ( nfac[i-1] ) - { - case 2 : - /*transform for factor of 2 (including rotation factor)*/ - - retVal = pre_fOf2Trans() ; - if ( retVal == 0 ) factorOf2Transform () ; - - break ; - - case 4 : - - kspnn = kspan ; - kspan = kspan >> 2 ; /*kspan /= 4 */ - - retVal = factorOf4Transform () ; - break ; - - case 3 : - - k = nfac[i-1] ; - kspnn = kspan ; - kspan = kspan / k ; - - factorOf3Transform ( ) ; - break ; - - case 5 : - - k = nfac[i-1] ; - kspnn = kspan ; - kspan = kspan / k ; - - factorOf5Transform ( ) ; - break ; - - default : - - k = nfac[i-1] ; - kspnn = kspan ; - kspan = kspan / k ; - - if ( nfac[i-1] != jf) preFOtherTransform ( ) ; - - factorOfOtherTransform ( ) ; - break ; - } - - - - if ( retVal == 42 ) - { - if ( i != m) retVal = mulByRotationFactor ( ) ; - else retVal = 1 ; - } - - if ( retVal == 1 ) return 1 ; /*goto permute */ - else return 0 ; /*goto factor_transform => once again*/ - -} - -/* permutation for square factor of n */ -static void permute_stage1 (void) -{ - - int retVal = 1 ; - - pre_sqFactor2NormlOrder () ; - - if ( n == ntot ) - /*permutation for single-variate transform (optional code)*/ - while ( retVal == 1) - { - single_sqFactor2NormlOrder () ; - retVal = post_sqFactor2NormlOrder () ; - } - else - /*permutation for multivariate transform*/ - while ( retVal == 1) retVal = multi_sqFactor2NormlOrder (); - -} - -static void permute_stage2 (void) -{ - kspnn = np[kt] ; - - /*permutation for square-free facotrs of n */ - nonSqFactor2NormOrder () ; - - /*determine the permutation cycles of length greater than 1*/ - detPermutCycles (); - - j = k3 + 1; - nt -= kspnn ; - i = nt - inc + 1 ; - while ( nt >= 0 ) - { - reorderMatrix ( ) ; - - j = k3 + 1 ; - nt -= kspnn ; - i = nt - inc + 1 ; - } -} - -/***************************************** -Sous-Sous-Fonctions -******************************************/ - - - - - -static int pre_fOf2Trans (void) -{ - kspan /= 2; - k1 = kspan + 2 ; - /*50*/ - do{ - do{ - k2 = kk + kspan ; - ak = a[k2-1] ; - bk = b[k2-1] ; - - a[k2-1] = a[kk-1] - ak; - b[k2-1] = b[kk-1] - bk; - - a[kk-1] = a[kk-1] + ak; - b[kk-1] = b[kk-1] + bk; - - kk = k2 + kspan ; - }while (kk <= nn); - - kk -= nn ; - }while (kk <= jc); - - - if ( kk > kspan ) return 1 ; /*goto350*/ - else return 0 ; /*goto60*/ - - -} - - - -static int factorOf2Transform (void) -{ - do /*60*/ {/*while ( kk <= jc*2 )*/ - c1 = 1 - cd ; - s1 = sd ; - mm = min( k1/2 , klim); - - do/* do 80 */ {/*while ( kk <= mm || ( kk > mm && kk < k2 ))*/ - do {/*while(kk > k2) */ - do { /*while ( kk < nt )*/ - k2 = kk + kspan; - - ak = a[kk-1] - a[k2-1]; - bk = b[kk-1] - b[k2-1]; - - a[kk-1] = a[kk-1] + a[k2-1]; - b[kk-1] = b[kk-1] + b[k2-1]; - - a[k2-1] = c1*ak - s1*bk; - b[k2-1] = s1*ak + c1*bk; - - kk = k2 + kspan; - }while ( kk < nt ); - - k2 = kk - nt; - c1 = -c1; - kk = k1 - k2; - - - }while (kk > k2); - - kk += jc; - - if ( kk <= mm ) /* 70 */ - { - ak = c1 - ( cd*c1+sd*s1) ; - s1 += (sd*c1-cd*s1) ; - /*c the following three statements compensate for truncation - c error. if rounded arithmetic is used, substitute - c c1=ak*/ - c1 = 0.5/(ak*ak+s1*s1) + 0.5 ; - s1 *= c1 ; - c1 *= ak ; - } - else { - if ( kk < k2 ) /*90*/ { - s1 = dr*rad*((double)(kk-1)/(double)jc); - c1 = cos(s1) ; - s1 = sin(s1) ; - mm = min(k1/2,mm+klim); - } - } - - } while ( kk <= mm || ( kk > mm && kk < k2 )); - - k1 += (inc+inc) ; - kk = (k1-kspan)/2 + jc; - - } while ( kk <= jc*2 ); - - - return 0 ; /*goto40*/ -} - - -/* this one is just an optimisation of the factor of 2 transform , we compute more things each turn */ - -static int factorOf4Transform (void) -{ - int return_value = 0 ; - - /*120*/ - do - { - c1 = 1 ; - s1 = 0 ; - - mm = min ( kspan , klim ) ; - - do - { - f4t_150 () ; - return_value = f4t_170 () ; - } while ( return_value == 0 ); - - kk += ( inc - kspan ) ; - } while ( kk <= jc ) ; - - - - if ( kspan == jc ) - return 1 ; /*goto350*/ - else - return 0 ;/*goto40*/ - - -} - -/*this function and the following are just here for conveniance , they just do fourier transformation for factor of 4 - but as the code was a bit long in factorof4transform , we've created two sub-functions */ - -static void f4t_150 (void) -{ - - do{ - k1 = kk + kspan ; - k2 = k1 + kspan ; - k3 = k2 + kspan ; - - akp = a[kk-1] + a[k2-1] ; - akm = a[kk-1] - a[k2-1] ; - - ajp = a[k1-1] + a[k3-1] ; - ajm = a[k1-1] - a[k3-1] ; - - a[kk-1] = akp + ajp ; - ajp = akp - ajp ; - - bkp = b[kk-1] + b[k2-1] ; - bkm = b[kk-1] - b[k2-1] ; - - bjp = b[k1-1] + b[k3-1] ; - bjm = b[k1-1] - b[k3-1] ; - - b[kk-1] = bkp + bjp ; - bjp = bkp - bjp ; - - akp = akm + bjm ; - akm = akm - bjm ; - - bkp = bkm - ajm ; - bkm = bkm + ajm ; - - if ( s1 == 0 )/*190*/ - { - a[k1-1] = akp ; - a[k2-1] = ajp ; - a[k3-1] = akm ; - - b[k1-1] = bkp ; - b[k2-1] = bjp ; - b[k3-1] = bkm ; - - - } - else /*160*/ - { - - a[k1-1] = akp*c1 - bkp*s1 ; - a[k2-1] = ajp*c2 - bjp*s2 ; - a[k3-1] = akm*c3 - bkm*s3 ; - - a[k1-1] = bkp*c1 + akp*s1 ; - a[k2-1] = bjp*c2 + ajp*s2 ; - a[k3-1] = bkm*c3 + akm*s3 ; - } - kk=k3+kspan; - }while ( kk <= nt ) ; - - -} - -static int f4t_170 (void) -{ - kk += ( jc - nt ) ; - - if ( kk <= mm ) - { - c2 = c1 - (cd*c1 + sd*s1); - s1 = s1 + (sd*c1 - cd*s1); - - /* - the following three statements compensate for truncation - error. if rounded arithmetic is used, substitute - c1=c2 - */ - - c1 = 0.5/(c2*c2+s1*s1) + 0.5 ; - s1 *= c1 ; - c1 *= c2 ; - - /*140*/ - - c2 = c1*c1 - s1*s1 ; - s2 = c1*s1*2 ; - c3 = c2*c1 - s2*s1 ; - s3 = c2*s1 + s2*c1 ; - - - return 0 ; - - } - else - { - if ( kk <= kspan ) - { - s1 = dr*rad * (kk-1)/jc ; - c1 = cos (s1) ; - s1 = sin (s1) ; - mm = min ( kspan , mm + klim ); - - /*140*/ - - c2 = c1*c1 - s1*s1 ; - s2 = c1*s1*2 ; - c3 = c2*c1 - s2*s1 ; - s3 = c2*s1 + s2*c1 ; - - return 0 ; - } - } - - return 1 ; -} - - - - -static void factorOf3Transform (void) -{ - do{ - do{ - k1 = kk + kspan ; - k2 = k1 + kspan ; - - ak = a[kk-1] ; - bk = b[kk-1] ; - - aj = a[k1-1] + a[k2-1] ; - bj = b[k1-1] + b[k2-1] ; - - a[kk-1] = ak + aj ; - b[kk-1] = bk + bj ; - - ak = -0.5*aj + ak ; - bk = -0.5*bj + bk ; - - aj = (a[k1-1] - a[k2-1])*s120 ; - bj = (b[k1-1] - b[k2-1])*s120 ; - - a[k1-1] = ak - bj ; - b[k1-1] = bk + aj ; - a[k2-1] = ak + bj ; - b[k2-1] = bk - aj ; - - kk = k2 + kspan ; - } while (kk < nn); - - kk -= nn ; - }while (kk <= kspan); - -} - -static void factorOf5Transform (void) -{ - c2 = c72*c72 - s72 *s72 ; - s2 = 2 * c72*s72; - - do{ - do{ - k1 = kk + kspan ; - k2 = k1 + kspan ; - k3 = k2 + kspan ; - k4 = k3 + kspan ; - - - - akp = a[k1-1] + a[k4-1] ; - akm = a[k1-1] - a[k4-1] ; - - bkp = b[k1-1] + b[k4-1] ; - bkm = b[k1-1] - b[k4-1] ; - - ajp = a[k2-1] + a[k3-1] ; - ajm = a[k2-1] - a[k3-1] ; - - bjp = b[k2-1] + b[k3-1] ; - bjm = b[k2-1] - b[k3-1] ; - - aa = a[kk-1] ; - bb = b[kk-1] ; - - a[kk-1] = aa + akp + ajp; - b[kk-1] = bb + bkp + bjp; - - ak = akp*c72 + ajp*c2 + aa ; - bk = bkp*c72 + bjp*c2 + bb ; - - aj = akm*s72 + ajm*s2 ; - bj = bkm*s72 + bjm*s2 ; - - a[k1-1] = ak - bj ; - a[k4-1] = ak + bj ; - b[k1-1] = bk + aj ; - b[k4-1] = bk - aj ; - - ak = akp*c2 + ajp*c72 + aa ; - bk = bkp*c2 + bjp*c72 + bb ; - - aj = akm*s2 - ajm*s72 ; - - bj = bkm*s2 - bjm*s72 ; - - a[k2-1] = ak - bj ; - a[k3-1] = ak + bj ; - b[k2-1] = bk + aj ; - b[k3-1] = bk - aj ; - - kk = k4 + kspan; - }while (kk < nn); - - kk -= nn ; - }while (kk <= kspan); -} - -/* this function is the general case of non factor of 2 factor , the factorof3transform and factorof5trandform are just -special case of this one */ - - -static void preFOtherTransform (void) -{ - - jf = k ; - s1 = (rad*8)/k ; - c1 = cos (s1) ; - s1 = sin (s1) ; - ck[jf-1] = 1 ; - sk[jf-1] = 0 ; - j = 1 ; - - do - { - ck[j-1] = ck[k-1] * c1 + sk[k-1]*s1 ; - sk[j-1] = ck[k-1] * s1 - sk[k-1]*c1 ; - - k -- ; - - ck[k-1] = ck[j-1] ; - sk[k-1] = - sk[j-1] ; - - j++ ; - - }while ( j < k ); - -} - -static void factorOfOtherTransform (void) -{ -int ktemp = 0 ; - -do - { - k1 = kk ; - k2 = kk + kspnn ; - - aa = a[kk-1] ; - bb = b[kk-1] ; - - ak = aa ; - bk = bb ; - j = 1 ; - k1 += kspan ; - - do - { - k2 -= kspan ; - j++ ; - - wt[j-1] = a[k1-1] + a[k2-1] ; - ak = wt[j-1] + ak ; - - bt[j-1] = b[k1-1] + b[k2-1] ; - bk = bt[j-1] + bk ; - j++ ; - - wt[j-1] = a[k1-1] - a[k2-1] ; - bt[j-1] = b[k1-1] - b[k2-1] ; - - k1 += kspan; - }while ( k1 < k2 ) ; - - a[kk-1] = ak ; - b[kk-1] = bk ; - - k1 = kk ; - k2 = kk + kspnn ; - j = 1 ; - - do - { - k1 += kspan ; - k2 -= kspan ; - jj = j ; - ak = aa ; - bk = bb ; - aj = 0 ; - bj = 0 ; - k = 1 ; - - do - { - k++ ; - ak += ( wt[k-1] * ck[jj-1] ) ; - bk += ( bt[k-1] * ck[jj-1] ) ; - - k++ ; - aj += (wt[k-1] * sk[jj-1]) ; - bj += (bt[k-1] * sk[jj-1]) ; - jj += j ; - - if ( jj > jf ) - jj -= jf ; - } while ( k < jf ) ; - - k = jf - j ; - a[k1-1] = ak - bj ; - b[k1-1] = bk + aj ; - a[k2-1] = ak + bj ; - b[k2-1] = bk - aj ; - - j++ ; - - }while ( j < k ) ; - - - - - - kk += kspnn ; - ktemp = kk ; - - if ( kk > nn ) - { - kk -= nn; - } - -}while ( ktemp <= nn || (kk <= kspan && ( ktemp > nn)) ); -} - - - - -static int mulByRotationFactor (void ) -{ - int ktemp = 0 ; - - if ( i != m ) - { - kk = jc + 1 ; - - /*300*/ - do - { - c2 = 1 - cd ; - s1 = sd ; - - mm = min ( kspan , klim ) ; - - /*320 */ - do - { - - c1 = c2 ; - s2 = s1 ; - kk += kspan ; - - do - { - - ak = a[kk-1] ; - a[kk-1] = c2*ak - s2*b[kk-1] ; - b[kk-1] = s2*ak + c2*b[kk-1] ; - - kk += kspnn ; - ktemp = kk ; - - - if ( kk > nt ) - { - ak = s1*s2 ; - s2 = s1*c2 + s2*c1 ; - c2 = c1*c2 - ak ; - kk += (kspan - nt ) ; - - - - } - - }while (ktemp <= nt || ( kk <= kspnn && ktemp > nt )) ; - - kk += ( jc - kspnn ); - - - if ( kk <= mm ) - { - - /* 310*/ - c2 = c1 - ( cd*c1 + sd*s1 ) ; - s1 += (sd*c1 - cd*s1 ) ; - - /* - the following three statements compensate for truncation - error. if rounded arithmetic is used, substitute - c1=c2 - */ - - c1 = 0.5/(c2*c2+s1*s1) + 0.5 ; - s1 *= c1 ; - c2 *= c1 ; - } - else - { - if ( kk <= kspan ) - { - - s1 = dr*rad * (kk-1)/jc ; - c2 = cos (s1) ; - s1 = sin (s1) ; - mm = min ( kspan , mm + klim ); - } - } - - }while ( kk <= mm ||( kk <= kspan && kk > mm ) ) ; - - kk += (jc + inc -kspan ); - - - }while ( kk <= jc+jc); - - - return 0 ; /* goto40 */ - } - - return 1 ; /* goto350*/ -} - - - - -static void pre_sqFactor2NormlOrder (void) -{ - - k = kt + kt + 1 ; - - if ( m < k ) - k -- ; - - j = 1 ; - np[k] = jc ; - - do - { - np[j] = np[j-1]/nfac[j-1] ; - np[k-1] = np[k]*nfac[j-1] ; - - j++ ; - k-- ; - }while ( j < k ) ; - - k3 = np[k] ; - kspan = np[1] ; - kk = jc + 1 ; - k2 = kspan + 1 ; - j = 1; - - -} - -static int post_sqFactor2NormlOrder (void) -{ - - do - { - do - { - k2 -= np[j-1] ; - j++ ; - k2 += np[j] ; - - } while ( k2 > np[j-1]); - - j = 1 ; - -/* 390 */ - do - { - if ( kk < k2 ) - { - - - return 1 ; - } - else - { - kk += inc ; - k2 += kspan ; - } - }while( k2 < ks ); - - }while ( kk < ks ) ; - - jc = k3 ; - - return 0; -} - - -/* appeler cetter fonction dans un do while valeur_retour != 1)*/ -static void single_sqFactor2NormlOrder (void) -{ - - - do - { - - ak = a[kk-1] ; - a[kk-1] = a[k2-1] ; - a[k2-1] = ak ; - - bk = b[kk-1] ; - b[kk-1] = b[k2-1] ; - b[k2-1] = bk ; - - kk += inc ; - k2 += kspan ; - } while ( k2 < ks ); - - /*380*/ -} - -/*idem que single_ */ -static int multi_sqFactor2NormlOrder (void) -{ - - - - k = kk + jc ; - - do /*410*/ - { - ak = a[kk-1] ; - a[kk-1] = a[k2-1] ; - a[k2-1] = ak ; - bk = b[kk-1] ; - b[kk-1] = b[k2-1] ; - b[k2-1] = bk ; - kk += inc ; - k2 += kspan ; - } while ( kk < k ); - - kk += (ks - jc ) ; - k2 += (ks - jc ) ; - - if ( kk < nt ) - return 1 ; - - k2 += ( kspan - nt ); - kk += ( jc - nt ); - - if ( k2 < ks ) - { - - return 1 ; - } - if( post_sqFactor2NormlOrder ( ) == 1 ) - { - - return 1 ; - } - jc = k3 ; - - return 0; - -} - - - -static void nonSqFactor2NormOrder (void) -{ - - j = m - kt ; - nfac[j] = 1 ; - - - - do - { - nfac[j-1] *= nfac[j] ; - - j-- ; - - - }while ( j != kt ) ; - - kt ++ ; - nn = nfac[kt-1] - 1; - - jj = 0 ; - j = 0; - - /*480*/ - - k2 = nfac[kt-1] ; - k = kt + 1 ; - kk = nfac[k-1] ; - j ++ ; - - while ( j <= nn ) - { - jj += kk ; - - - while ( jj >= k2 ) - { - jj -= k2 ; - k2 = kk ; - k++ ; - kk = nfac[k-1] ; - - jj += kk ; - - - } - - np[j-1] = jj ; - k2 = nfac[kt-1] ; - k = kt + 1 ; - kk = nfac[k-1] ; - j ++ ; - - } - - j = 0 ; - - return ; -} - -/* here we determine how many permutation cycles we need to do */ -static void detPermutCycles (void) -{ - - do - { - do - { - j++ ; - - kk = np[j-1] ; - }while ( kk < 0 ) ; - - - if ( kk != j ) - { - do - { - - k = kk ; - kk = np[k-1] ; - np[k-1] = -kk ; - }while ( kk != j ) ; - k3 = kk ; - } - else - np[j-1] = -j ; - }while ( j != nn ); - - maxf *= inc ; - - return ; -} - -static void reorderMatrix (void) -{ -do - { - do - { - j-- ; - - }while (np[j-1] < 0 ) ; - - jj = jc ; - - /*520*/ - do - { - kspan = jj ; - - if ( jj > maxf ) - kspan = maxf ; - - jj -= kspan ; - k = np [j-1]; - kk = jc*k + i + jj ; - k1 = kk + kspan ; - k2 = 0 ; - - do /*530*/ - { - k2 ++ ; - wt[k2-1] = a[k1-1] ; - bt[k2-1] = b[k1-1] ; - k1 -= inc ; - - }while ( k1 != kk ); - - do - { - k1 = kk + kspan ; - k2 = k1 - jc * (k + np[k-1]); - k = -np[k-1]; - - - do - { - a[k1-1] = a[k2-1] ; - b[k1-1] = b[k2-1] ; - - k1 -= inc ; - k2 -= inc ; - - }while ( k1 != kk ) ; - - kk = k2 ; - - }while ( k != j ); - - k1 = kk +kspan ; - k2 = 0 ; - /*560*/ - - do - { - k2 ++ ; - a[k1-1] = wt[k2-1] ; - b[k1-1] = bt[k2-1] ; - k1 -= inc ; - - - }while ( k1 != kk ) ; - - } while ( jj != 0 ) ; -}while ( j != 1 ) ; - - return ; -} - - -- cgit