@@ -484,6 +484,7 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
484
484
485
485
sotDEBUGF (5 , " --- Time %d -------------------" , iterTime );
486
486
unsigned int iterTask = 0 ;
487
+ const Matrix* PrevProj = NULL ;
487
488
for ( StackType::iterator iter = stack.begin (); iter!=stack.end ();++iter )
488
489
{
489
490
sotDEBUGF (5 ," Rank %d." ,iterTask);
@@ -504,17 +505,16 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
504
505
task.memoryInternal = mem;
505
506
}
506
507
507
- dynamicgraph:: Matrix &Jp = mem->Jp ;
508
- dynamicgraph:: Matrix &V = mem->V ;
509
- dynamicgraph:: Matrix &JK = mem->JK ;
510
- dynamicgraph:: Matrix &Jt = mem->Jt ;
508
+ Matrix &Jp = mem->Jp ;
509
+ Matrix &JK = mem->JK ;
510
+ Matrix &Jt = mem->Jt ;
511
+ Matrix &Proj = mem->Proj ;
511
512
MemoryTaskSOT::SVD_t& svd = mem->svd ;
512
513
513
514
taskVectorToMlVector (task.taskSOUT (iterTime), mem->err );
514
515
const dynamicgraph::Vector &err = mem->err ;
515
516
516
517
Jp.resize ( mJ ,nJ );
517
- V.resize ( mJ ,mJ );
518
518
Jt.resize ( nJ,mJ );
519
519
JK.resize ( nJ,mJ );
520
520
@@ -551,7 +551,7 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
551
551
/* **/ sotCOUNTER (2 ,3 ); // compute JK
552
552
553
553
/* --- COMPUTE Jt --- */
554
- if ( 0 <iterTask ) Jt.noalias () = JK*Proj ; else { Jt = JK; }
554
+ if ( 0 <iterTask ) Jt.noalias () = JK*(*PrevProj) ; else { Jt = JK; }
555
555
/* **/ sotCOUNTER (3 ,4 ); // compute Jt
556
556
557
557
/* --- COMPUTE S --- */
@@ -561,9 +561,8 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
561
561
/* --- PINV --- */
562
562
svd.compute (Jt);
563
563
Eigen::dampedInverse (svd, Jp, th);
564
- V.noalias () = svd.matrixV ();
565
564
/* **/ sotCOUNTER (5 ,6 ); // PINV
566
- sotDEBUG (2 ) << " V after dampedInverse." << V <<endl;
565
+ sotDEBUG (20 ) << " V after dampedInverse." << svd. matrixV () <<endl;
567
566
/* --- RANK --- */
568
567
{
569
568
rankJ = 0 ;
@@ -581,7 +580,7 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
581
580
sotDEBUG (45 ) << " JJp" <<iterTask<<" = " << JK*Jp <<endl;
582
581
// sotDEBUG(45) << "U"<<iterTask<<" = "<< U<<endl;
583
582
sotDEBUG (45 ) << " S" <<iterTask<<" = " << svd.singularValues ()<<endl;
584
- sotDEBUG (45 ) << " V" <<iterTask<<" = " << V <<endl;
583
+ sotDEBUG (45 ) << " V" <<iterTask<<" = " << svd. matrixV () <<endl;
585
584
sotDEBUG (45 ) << " U" <<iterTask<<" = " << svd.matrixU ()<<endl;
586
585
587
586
mem->jacobianInvSINOUT = Jp;
@@ -590,7 +589,7 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
590
589
mem->jacobianConstrainedSINOUT .setTime ( iterTime );
591
590
mem->jacobianProjectedSINOUT = Jt;
592
591
mem->jacobianProjectedSINOUT .setTime ( iterTime );
593
- mem->singularBaseImageSINOUT = V ;
592
+ mem->singularBaseImageSINOUT = svd. matrixV (). leftCols (rankJ) ;
594
593
mem->singularBaseImageSINOUT .setTime ( iterTime );
595
594
mem->rankSINOUT = rankJ;
596
595
mem->rankSINOUT .setTime ( iterTime );
@@ -600,27 +599,23 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
600
599
sotDEBUG (2 )<<" Inverse not recomputed." <<endl;
601
600
rankJ = mem->rankSINOUT .accessCopy ();
602
601
Jp = mem->jacobianInvSINOUT .accessCopy ();
603
- V = mem->singularBaseImageSINOUT .accessCopy ();
604
602
JK = mem->jacobianConstrainedSINOUT .accessCopy ();
605
603
Jt = mem->jacobianProjectedSINOUT .accessCopy ();
606
604
}
607
605
/* **/ sotCOUNTER (6 ,7 ); // TRACE
608
606
609
- if (rankJ == 0 ) break ;
610
-
611
607
/* --- COMPUTE QDOT AND P --- */
612
608
/* DEBUG: normally, the first iter (ie the test below)
613
609
* is the same than the other, starting with control_0 = q0SIN. */
614
610
if ( iterTask==0 ) control.noalias () += Jp*err;
615
- else control += Proj * (Jp*(err - JK*control));
611
+ else control += *PrevProj * (Jp*(err - JK*control));
616
612
/* **/ sotCOUNTER (7 ,8 ); // QDOT
617
613
618
614
/* --- OPTIMAL FORM: To debug. --- */
619
615
if ( 0 ==iterTask ) {
620
- Proj.resize ( mJ , mJ ); Proj. setIdentity ( );
616
+ Proj.noalias () = svd. matrixV (). rightCols (svd. matrixV (). cols ()-rankJ );
621
617
} else {
622
- // Proj.noalias() -= svd.matrixV().leftCols(rankJ) * svd.matrixV().leftCols(rankJ).adjoint();
623
- Proj = Proj * svd.matrixV ().rightCols (svd.matrixV ().cols ()-rankJ);
618
+ Proj.noalias () = *PrevProj * svd.matrixV ().rightCols (svd.matrixV ().cols ()-rankJ);
624
619
}
625
620
626
621
/* --- OLIVIER START --- */
@@ -629,16 +624,17 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
629
624
sotDEBUG (2 ) << " Proj non optimal (rankJ= " <<rankJ
630
625
<< " , iterTask =" << iterTask
631
626
<< " )" ;
632
- sotDEBUG (2 ) << " V = " << V ;
633
- sotDEBUG (2 ) << " Jt = " << Jt;
634
- sotDEBUG (2 ) << " JpxJt = " << Jp*Jt;
627
+ sotDEBUG (20 ) << " V = " << svd. matrixV () ;
628
+ sotDEBUG (20 ) << " Jt = " << Jt;
629
+ sotDEBUG (20 ) << " JpxJt = " << Jp*Jt;
635
630
sotDEBUG (25 ) << " Proj-Jp*Jt" <<iterTask<<" = " << (Proj-Jp*Jt) <<endl;
636
631
637
632
/* --- OLIVIER END --- */
638
633
639
634
sotDEBUG (15 ) << " q" <<iterTask<<" = " <<control<<std::endl;
640
635
sotDEBUG (25 ) << " P" <<iterTask<<" = " << Proj <<endl;
641
636
iterTask++;
637
+ PrevProj = &Proj;
642
638
643
639
sotPRINTCOUNTER (1 ); sotPRINTCOUNTER (2 ); sotPRINTCOUNTER (3 );
644
640
sotPRINTCOUNTER (4 ); sotPRINTCOUNTER (5 ); sotPRINTCOUNTER (6 );
@@ -695,7 +691,10 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
695
691
svd.compute (Jt);
696
692
// TODO the two next lines could be replaced by
697
693
Eigen::dampedInverse ( svd, Jp,th );
698
- PJp.noalias () = Proj*Jp;
694
+ if (PrevProj != NULL )
695
+ PJp.noalias () = (*PrevProj)*Jp;
696
+ else
697
+ PJp.noalias () = Jp;
699
698
700
699
/* --- COMPUTE ERR --- */
701
700
const dynamicgraph::Vector& Herr ( err );
@@ -705,7 +704,7 @@ computeControlLaw( dynamicgraph::Vector& control,const int& iterTime )
705
704
706
705
/* --- TRACE --- */
707
706
sotDEBUG (45 ) << " Pgrad = " << (PJp*Herr) <<endl;
708
- sotDEBUG (45 ) << " P = " << Proj <<endl;
707
+ if (PrevProj != NULL ) { sotDEBUG (45 ) << " P = " << *PrevProj <<endl; }
709
708
sotDEBUG (45 ) << " Jp = " << Jp <<endl;
710
709
sotDEBUG (45 ) << " PJp = " << PJp <<endl;
711
710
}
0 commit comments