1 /*
2    CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
3    Copyright (C) 2013-2018 Sebastian Wouters
4 
5    This program is free software; you can redistribute it and/or modify
6    it under the terms of the GNU General Public License as published by
7    the Free Software Foundation; either version 2 of the License, or
8    (at your option) any later version.
9 
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU General Public License for more details.
14 
15    You should have received a copy of the GNU General Public License along
16    with this program; if not, write to the Free Software Foundation, Inc.,
17    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 
20 #include <iostream>
21 #include <stdlib.h>
22 #include <sstream>
23 #include <string>
24 #include <algorithm>
25 #include <math.h>
26 #include <assert.h>
27 #include <sys/time.h>
28 #include <unistd.h>
29 
30 #include "DMRG.h"
31 #include "Lapack.h"
32 #include "Heff.h"
33 #include "MPIchemps2.h"
34 #include "Excitation.h"
35 
36 using std::cout;
37 using std::endl;
38 using std::min;
39 using std::max;
40 
Symm4RDM(double * output,const int Y,const int Z,const bool last_case)41 void CheMPS2::DMRG::Symm4RDM( double * output, const int Y, const int Z, const bool last_case ){
42 
43    struct timeval start, end;
44    gettimeofday( &start, NULL );
45 
46    assert( the3DM != NULL );
47 
48    symm_4rdm_helper( output, Y, Z, 1.0, 1.0, false, 0.5 ); // output = 0.5 *   3rdm[ ( 1 + E_{YZ} + E_{ZY} ) | 0 > ]
49    symm_4rdm_helper( output, Y, Z, 1.0, 0.0, true, -0.5 ); // output = 0.5 * ( 3rdm[ ( 1 + E_{YZ} + E_{ZY} ) | 0 > ] - 3rdm[ E_{YZ} + E_{ZY} | 0 > ] )
50 
51    for ( int r = 0; r < L; r++ ){
52       for ( int q = 0; q < L; q++ ){
53          for ( int p = 0; p < L; p++ ){
54             for ( int k = 0; k < L; k++ ){
55                for ( int j = 0; j < L; j++ ){
56                   for ( int i = 0; i < L; i++ ){
57                      output[ i + L * ( j + L * ( k + L * ( p + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( i, j, k, p, q, r );
58                   }
59                   output[ Y + L * ( j + L * ( k + L * ( p + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( Z, j, k, p, q, r );
60                   output[ Z + L * ( j + L * ( k + L * ( p + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( Y, j, k, p, q, r );
61                   output[ j + L * ( Y + L * ( k + L * ( p + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( j, Z, k, p, q, r );
62                   output[ j + L * ( Z + L * ( k + L * ( p + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( j, Y, k, p, q, r );
63                   output[ j + L * ( k + L * ( Y + L * ( p + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( j, k, Z, p, q, r );
64                   output[ j + L * ( k + L * ( Z + L * ( p + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( j, k, Y, p, q, r );
65                   output[ j + L * ( k + L * ( p + L * ( Y + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( j, k, p, Z, q, r );
66                   output[ j + L * ( k + L * ( p + L * ( Z + L * ( q + L * r )))) ] -= 0.5 * the3DM->get_ham_index( j, k, p, Y, q, r );
67                   output[ j + L * ( k + L * ( p + L * ( q + L * ( Y + L * r )))) ] -= 0.5 * the3DM->get_ham_index( k, j, p, Z, q, r );
68                   output[ j + L * ( k + L * ( p + L * ( q + L * ( Z + L * r )))) ] -= 0.5 * the3DM->get_ham_index( k, j, p, Y, q, r );
69                   output[ j + L * ( k + L * ( p + L * ( q + L * ( r + L * Y )))) ] -= 0.5 * the3DM->get_ham_index( p, j, k, Z, q, r );
70                   output[ j + L * ( k + L * ( p + L * ( q + L * ( r + L * Z )))) ] -= 0.5 * the3DM->get_ham_index( p, j, k, Y, q, r );
71                }
72             }
73          }
74       }
75    }
76 
77    if ( last_case ){ PreSolve(); } // Need to set up the renormalized operators again to continue sweeping
78 
79    gettimeofday( &end, NULL );
80    const double elapsed = ( end.tv_sec - start.tv_sec ) + 1e-6 * ( end.tv_usec - start.tv_usec );
81    #ifdef CHEMPS2_MPI_COMPILATION
82    if ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER )
83    #endif
84    { cout << "CheMPS2::DMRG::Symm4RDM( " << Y << " , " << Z << " ) : Elapsed wall time = " << elapsed << " seconds." << endl; }
85 
86 }
87 
symm_4rdm_helper(double * output,const int ham_orb1,const int ham_orb2,const double alpha,const double beta,const bool add,const double factor)88 void CheMPS2::DMRG::symm_4rdm_helper( double * output, const int ham_orb1, const int ham_orb2, const double alpha, const double beta, const bool add, const double factor ){
89 
90    // Figure out the DMRG orbitals, in order
91    assert( ham_orb1 >= 0 );
92    assert( ham_orb2 >= 0 );
93    assert( ham_orb1 <  L );
94    assert( ham_orb2 <  L );
95    const int  first_dmrg_orb = (( Prob->gReorder() ) ? Prob->gf1( ham_orb1 ) : ham_orb1 );
96    const int second_dmrg_orb = (( Prob->gReorder() ) ? Prob->gf1( ham_orb2 ) : ham_orb2 );
97    const int dmrg_orb1 = (( first_dmrg_orb <= second_dmrg_orb ) ?  first_dmrg_orb : second_dmrg_orb );
98    const int dmrg_orb2 = (( first_dmrg_orb <= second_dmrg_orb ) ? second_dmrg_orb :  first_dmrg_orb );
99 
100    // Make a back-up of the entirely left-normalized MPS and the corresponding bookkeeper
101    SyBookkeeper * oldBK = denBK;
102    if ( dmrg_orb1 != dmrg_orb2 ){ denBK = new SyBookkeeper( *oldBK ); }
103    TensorT ** backup_mps = new TensorT * [ L ];
104    for ( int orbital = 0; orbital < L; orbital++ ){
105       backup_mps[ orbital ] = MPS[ orbital ];
106       MPS[ orbital ] = new TensorT( orbital, denBK ); // denBK is now a DIFFERENT pointer than backup_mps[ orbital ]->gBK()
107       int totalsize = MPS[ orbital ]->gKappa2index( MPS[ orbital ]->gNKappa() );
108       int inc1 = 1;
109       dcopy_( &totalsize, backup_mps[ orbital ]->gStorage(), &inc1, MPS[ orbital ]->gStorage(), &inc1 );
110    }
111    deleteAllBoundaryOperators();
112 
113    // Change the gauge so that the non-orthonormal MPS tensor is on site dmrg_orb2
114    for ( int siteindex = L - 1; siteindex > dmrg_orb2; siteindex-- ){
115       right_normalize( MPS[ siteindex - 1 ], MPS[ siteindex ] );
116       updateMovingLeftSafeFirstTime( siteindex - 1 );
117    }
118 
119    // Solve
120    solve_fock( dmrg_orb1, dmrg_orb2, alpha, beta );
121 
122    // Further right normalize the wavefunction except for the first MPS tensor ( contains the norm )
123    for ( int siteindex = dmrg_orb2; siteindex > 0; siteindex-- ){
124       right_normalize( MPS[ siteindex - 1 ], MPS[ siteindex ] );
125       updateMovingLeftSafeFirstTime( siteindex - 1 );
126    }
127 
128    ThreeDM * helper3rdm = new ThreeDM( denBK, Prob, false );
129    tensor_3rdm_a_J0_doublet = new Tensor3RDM****[ L - 1 ];
130    tensor_3rdm_a_J1_doublet = new Tensor3RDM****[ L - 1 ];
131    tensor_3rdm_a_J1_quartet = new Tensor3RDM****[ L - 1 ];
132    tensor_3rdm_b_J0_doublet = new Tensor3RDM****[ L - 1 ];
133    tensor_3rdm_b_J1_doublet = new Tensor3RDM****[ L - 1 ];
134    tensor_3rdm_b_J1_quartet = new Tensor3RDM****[ L - 1 ];
135    tensor_3rdm_c_J0_doublet = new Tensor3RDM****[ L - 1 ];
136    tensor_3rdm_c_J1_doublet = new Tensor3RDM****[ L - 1 ];
137    tensor_3rdm_c_J1_quartet = new Tensor3RDM****[ L - 1 ];
138    tensor_3rdm_d_J0_doublet = new Tensor3RDM****[ L - 1 ];
139    tensor_3rdm_d_J1_doublet = new Tensor3RDM****[ L - 1 ];
140    tensor_3rdm_d_J1_quartet = new Tensor3RDM****[ L - 1 ];
141 
142    // Leftmost contribution to the helper3rdm
143    helper3rdm->fill_site( MPS[ 0 ], Ltensors, F0tensors, F1tensors, S0tensors, S1tensors, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL );
144 
145    // Other contributions to the helper3rdm
146    for ( int siteindex = 1; siteindex < L; siteindex++ ){
147 
148       /* Change the MPS gauge */
149       left_normalize( MPS[ siteindex - 1 ], MPS[ siteindex ] );
150 
151       /* Update the required renormalized operators */
152       update_safe_3rdm_operators( siteindex );
153       updateMovingRightSafe2DM( siteindex - 1 );
154 
155       /* Current contribution to helper3rdm */
156       helper3rdm->fill_site( MPS[ siteindex ], Ltensors, F0tensors, F1tensors, S0tensors, S1tensors,
157                              tensor_3rdm_a_J0_doublet[ siteindex - 1 ], tensor_3rdm_a_J1_doublet[ siteindex - 1 ], tensor_3rdm_a_J1_quartet[ siteindex - 1 ],
158                              tensor_3rdm_b_J0_doublet[ siteindex - 1 ], tensor_3rdm_b_J1_doublet[ siteindex - 1 ], tensor_3rdm_b_J1_quartet[ siteindex - 1 ],
159                              tensor_3rdm_c_J0_doublet[ siteindex - 1 ], tensor_3rdm_c_J1_doublet[ siteindex - 1 ], tensor_3rdm_c_J1_quartet[ siteindex - 1 ],
160                              tensor_3rdm_d_J0_doublet[ siteindex - 1 ], tensor_3rdm_d_J1_doublet[ siteindex - 1 ], tensor_3rdm_d_J1_quartet[ siteindex - 1 ] );
161 
162    }
163 
164    // Collect all data
165    #ifdef CHEMPS2_MPI_COMPILATION
166    helper3rdm->mpi_allreduce();
167    #endif
168 
169    // Copy the contributions
170    helper3rdm->correct_higher_multiplicities();
171    delete_3rdm_operators( L - 1 );
172    delete [] tensor_3rdm_a_J0_doublet;
173    delete [] tensor_3rdm_a_J1_doublet;
174    delete [] tensor_3rdm_a_J1_quartet;
175    delete [] tensor_3rdm_b_J0_doublet;
176    delete [] tensor_3rdm_b_J1_doublet;
177    delete [] tensor_3rdm_b_J1_quartet;
178    delete [] tensor_3rdm_c_J0_doublet;
179    delete [] tensor_3rdm_c_J1_doublet;
180    delete [] tensor_3rdm_c_J1_quartet;
181    delete [] tensor_3rdm_d_J0_doublet;
182    delete [] tensor_3rdm_d_J1_doublet;
183    delete [] tensor_3rdm_d_J1_quartet;
184    helper3rdm->fill_ham_index( factor, add, output, 0, L );
185 
186    // Throw out the changed MPS and place back the original left-normalized MPS
187    for ( int orbital = 0; orbital < L; orbital++ ){
188       delete MPS[ orbital ];
189       MPS[ orbital ] = backup_mps[ orbital ];
190    }
191    delete [] backup_mps;
192    if ( dmrg_orb1 != dmrg_orb2 ){
193       delete denBK;
194       denBK = oldBK;
195    }
196    delete helper3rdm;
197    deleteAllBoundaryOperators();
198 
199 }
200 
solve_fock(const int dmrg_orb1,const int dmrg_orb2,const double alpha,const double beta)201 void CheMPS2::DMRG::solve_fock( const int dmrg_orb1, const int dmrg_orb2, const double alpha, const double beta ){
202 
203    #ifdef CHEMPS2_MPI_COMPILATION
204       const bool am_i_master = ( MPIchemps2::mpi_rank() == MPI_CHEMPS2_MASTER );
205    #else
206       const bool am_i_master = true;
207    #endif
208 
209    if ( dmrg_orb1 == dmrg_orb2 ){
210       MPS[ dmrg_orb1 ]->number_operator( 2 * alpha, beta ); // alpha * ( E_zz + E_zz ) + beta * 1
211       return;
212    }
213 
214    double sweep_inproduct = 0.0;
215 
216    if ( dmrg_orb1 + 1 == dmrg_orb2 ){
217       Sobject * newS = new Sobject( dmrg_orb1, denBK );
218       if ( am_i_master ){
219          Sobject * oldS = new Sobject( dmrg_orb1, denBK );
220          oldS->Join( MPS[ dmrg_orb1 ], MPS[ dmrg_orb2 ] );
221          sweep_inproduct = Excitation::matvec( denBK, denBK, dmrg_orb1, dmrg_orb2, alpha, alpha, beta, newS, oldS, NULL, NULL, NULL );
222          delete oldS;
223       }
224       // MPI_CHEMPS2_MASTER decomposes newS. Each MPI process returns the correct discarded_weight. Each MPI process has the new MPS tensors set.
225       const double discarded_weight = newS->Split( MPS[ dmrg_orb1 ], MPS[ dmrg_orb2 ], OptScheme->get_D( OptScheme->get_number() - 1 ), true, true );
226       delete newS;
227    }
228 
229    if ( dmrg_orb1 + 1 < dmrg_orb2 ){
230 
231       SyBookkeeper * newBK = denBK;
232       denBK = new SyBookkeeper( *newBK );
233       newBK->restart( dmrg_orb1 + 1, dmrg_orb2, OptScheme->get_D( OptScheme->get_number() - 1 ) );
234       TensorT ** old_mps = new TensorT * [ L ];
235       for ( int orbital = dmrg_orb1; orbital <= dmrg_orb2; orbital++ ){
236          old_mps[ orbital ] = MPS[ orbital ];
237          old_mps[ orbital ]->sBK( denBK );
238          MPS[ orbital ] = new TensorT( orbital, newBK );
239          MPS[ orbital ]->random();
240          left_normalize( MPS[ orbital ], NULL ); // MPI_CHEMPS2_MASTER broadcasts MPS[ orbital ] ( left-normalized ).
241       }
242 
243       TensorO ** overlaps = NULL;
244       TensorL ** regular  = NULL;
245       TensorL ** trans    = NULL;
246       if ( am_i_master ){
247          overlaps = new TensorO*[ L - 1 ];
248          regular  = new TensorL*[ L - 1 ];
249          trans    = new TensorL*[ L - 1 ];
250          for ( int cnt = 0; cnt < L - 1; cnt++ ){
251             overlaps[ cnt ] = NULL;
252              regular[ cnt ] = NULL;
253                trans[ cnt ] = NULL;
254          }
255          for ( int orbital = dmrg_orb1; orbital < dmrg_orb2 - 1; orbital++ ){
256             solve_fock_update_helper( orbital, dmrg_orb1, dmrg_orb2, true, MPS, old_mps, newBK, denBK, overlaps, regular, trans );
257          }
258       }
259 
260       // Sweeps
261       bool change = false;
262       for ( int instruction = 0; instruction < OptScheme->get_number(); instruction++ ){
263          int num_iterations = 0;
264          double previous_inproduct = sweep_inproduct + 10 * OptScheme->get_energy_conv( instruction );
265          while (( fabs( sweep_inproduct - previous_inproduct ) > OptScheme->get_energy_conv( instruction ) ) && ( num_iterations < OptScheme->get_max_sweeps( instruction ) )){
266             {
267                const double noise_level = fabs( OptScheme->get_noise_prefactor( instruction ) ) * MaxDiscWeightLastSweep;
268                MaxDiscWeightLastSweep   = 0.0;
269                for ( int index = dmrg_orb2 - 1; index > dmrg_orb1; index-- ){
270                   Sobject * newS = new Sobject( index, newBK );
271                   if ( am_i_master ){
272                      Sobject * oldS = new Sobject( index, denBK );
273                      oldS->Join( old_mps[ index ], old_mps[ index + 1 ] );
274                      sweep_inproduct = Excitation::matvec( newBK, denBK, dmrg_orb1, dmrg_orb2, alpha, alpha, beta, newS, oldS, overlaps, regular, trans );
275                      delete oldS;
276                      if ( noise_level > 0.0 ){ newS->addNoise( noise_level ); }
277                   }
278                   // MPI_CHEMPS2_MASTER decomposes newS. Each MPI process returns the correct discarded_weight. Each MPI process has the new MPS tensors set.
279                   const double discarded_weight = newS->Split( MPS[ index ], MPS[ index + 1 ], OptScheme->get_D( instruction ), false, change );
280                   if ( discarded_weight > MaxDiscWeightLastSweep ){ MaxDiscWeightLastSweep = discarded_weight; }
281                   delete newS;
282                   if ( am_i_master ){
283                      solve_fock_update_helper( index, dmrg_orb1, dmrg_orb2, false, MPS, old_mps, newBK, denBK, overlaps, regular, trans );
284                   }
285                }
286             }
287             change = true;
288             {
289                const double noise_level = fabs( OptScheme->get_noise_prefactor( instruction ) ) * MaxDiscWeightLastSweep;
290                MaxDiscWeightLastSweep   = 0.0;
291                for ( int index = dmrg_orb1; index < dmrg_orb2 - 1; index++ ){
292                   Sobject * newS = new Sobject( index, newBK );
293                   if ( am_i_master ){
294                      Sobject * oldS = new Sobject( index, denBK );
295                      oldS->Join( old_mps[ index ], old_mps[ index + 1 ] );
296                      sweep_inproduct = Excitation::matvec( newBK, denBK, dmrg_orb1, dmrg_orb2, alpha, alpha, beta, newS, oldS, overlaps, regular, trans );
297                      delete oldS;
298                      if ( noise_level > 0.0 ){ newS->addNoise( noise_level ); }
299                   }
300                   // MPI_CHEMPS2_MASTER decomposes newS. Each MPI process returns the correct discarded_weight. Each MPI process has the new MPS tensors set.
301                   const double discarded_weight = newS->Split( MPS[ index ], MPS[ index + 1 ], OptScheme->get_D( instruction ), true, change );
302                   if ( discarded_weight > MaxDiscWeightLastSweep ){ MaxDiscWeightLastSweep = discarded_weight; }
303                   delete newS;
304                   if ( am_i_master ){
305                      solve_fock_update_helper( index, dmrg_orb1, dmrg_orb2, true, MPS, old_mps, newBK, denBK, overlaps, regular, trans );
306                   }
307                }
308             }
309             #ifdef CHEMPS2_MPI_COMPILATION
310             CheMPS2::MPIchemps2::broadcast_array_double( &sweep_inproduct, 1, MPI_CHEMPS2_MASTER );
311             #endif
312             previous_inproduct = sweep_inproduct;
313             num_iterations++;
314          }
315       }
316 
317       if ( am_i_master ){
318          for ( int index = 0; index < L - 1; index++ ){
319             if ( overlaps[ index ] != NULL ){ delete overlaps[ index ]; }
320             if (  regular[ index ] != NULL ){ delete  regular[ index ]; }
321             if (    trans[ index ] != NULL ){ delete    trans[ index ]; }
322          }
323          delete [] overlaps;
324          delete [] regular;
325          delete [] trans;
326       }
327 
328       for ( int orbital = dmrg_orb1; orbital <= dmrg_orb2; orbital++ ){ delete old_mps[ orbital ]; }
329       delete [] old_mps;
330       delete denBK;
331       denBK = newBK;
332 
333    }
334 
335    if ( am_i_master ){
336       const double rdm_inproduct = 2 * alpha * the2DM->get1RDM_DMRG( dmrg_orb1, dmrg_orb2 ) + beta;
337       cout << "DMRG::solve_fock : Accuracy = " << fabs( sweep_inproduct / ( Prob->gTwoS() + 1 ) - rdm_inproduct ) << endl;
338    }
339 
340 }
341 
solve_fock_update_helper(const int index,const int dmrg_orb1,const int dmrg_orb2,const bool moving_right,TensorT ** new_mps,TensorT ** old_mps,SyBookkeeper * new_bk,SyBookkeeper * old_bk,TensorO ** overlaps,TensorL ** regular,TensorL ** trans)342 void CheMPS2::DMRG::solve_fock_update_helper( const int index, const int dmrg_orb1, const int dmrg_orb2, const bool moving_right, TensorT ** new_mps, TensorT ** old_mps, SyBookkeeper * new_bk, SyBookkeeper * old_bk, TensorO ** overlaps, TensorL ** regular, TensorL ** trans ){
343 
344    if ( overlaps[ index ] != NULL ){ delete overlaps[ index ]; }
345    if (  regular[ index ] != NULL ){ delete  regular[ index ]; }
346    if (    trans[ index ] != NULL ){ delete    trans[ index ]; }
347 
348    const int Idiff = new_bk->gIrrep( dmrg_orb1 );
349    assert( Idiff == new_bk->gIrrep( dmrg_orb2 ) );
350    overlaps[ index ] = new TensorO( index + 1,        moving_right, new_bk, old_bk );
351     regular[ index ] = new TensorL( index + 1, Idiff, moving_right, new_bk, old_bk );
352       trans[ index ] = new TensorL( index + 1, Idiff, moving_right, old_bk, new_bk );
353 
354    if ( moving_right ){
355       if ( index == dmrg_orb1 ){
356          overlaps[ index ]->create( new_mps[ index ], old_mps[ index ] );
357           regular[ index ]->create( new_mps[ index ], old_mps[ index ], NULL, NULL );
358             trans[ index ]->create( old_mps[ index ], new_mps[ index ], NULL, NULL );
359       } else {
360          const int dimL = std::max( new_bk->gMaxDimAtBound( index     ), old_bk->gMaxDimAtBound( index     ) );
361          const int dimR = std::max( new_bk->gMaxDimAtBound( index + 1 ), old_bk->gMaxDimAtBound( index + 1 ) );
362          double * workmem = new double[ dimL * dimR ];
363          overlaps[ index ]->update( overlaps[ index - 1 ], new_mps[ index ], old_mps[ index ], workmem );
364           regular[ index ]->update(  regular[ index - 1 ], new_mps[ index ], old_mps[ index ], workmem );
365             trans[ index ]->update(    trans[ index - 1 ], old_mps[ index ], new_mps[ index ], workmem );
366          delete [] workmem;
367       }
368    } else {
369       if ( index + 1 == dmrg_orb2 ){
370          overlaps[ index ]->create( new_mps[ index + 1 ], old_mps[ index + 1 ] );
371           regular[ index ]->create( new_mps[ index + 1 ], old_mps[ index + 1 ], NULL, NULL );
372             trans[ index ]->create( old_mps[ index + 1 ], new_mps[ index + 1 ], NULL, NULL );
373       } else {
374          const int dimL = std::max( new_bk->gMaxDimAtBound( index + 1 ), old_bk->gMaxDimAtBound( index + 1 ) );
375          const int dimR = std::max( new_bk->gMaxDimAtBound( index + 2 ), old_bk->gMaxDimAtBound( index + 2 ) );
376          double * workmem = new double[ dimL * dimR ];
377          overlaps[ index ]->update( overlaps[ index + 1 ], new_mps[ index + 1 ], old_mps[ index + 1 ], workmem );
378           regular[ index ]->update(  regular[ index + 1 ], new_mps[ index + 1 ], old_mps[ index + 1 ], workmem );
379             trans[ index ]->update(    trans[ index + 1 ], old_mps[ index + 1 ], new_mps[ index + 1 ], workmem );
380          delete [] workmem;
381       }
382    }
383 
384 }
385 
386 
387