diff --git a/devel/cc/curly_Fock.irp.f b/devel/cc/curly_Fock.irp.f
index aecff09..0f22e3a 100644
--- a/devel/cc/curly_Fock.irp.f
+++ b/devel/cc/curly_Fock.irp.f
@@ -6,36 +6,41 @@ BEGIN_PROVIDER [ double precision, c_spin_fock_matrix_mo_oo, (spin_occ_num,spin_
   integer                        :: i,j,m,n
   integer                        :: a,b,e,f
   
-  do m=1,spin_occ_num
-    do i=1,spin_occ_num
-      
+  do i=1,spin_occ_num
+    do m=1,spin_occ_num
       if (m /= i) then
         c_spin_fock_matrix_mo_oo(m,i) = spin_fock_matrix_mo_oo(m,i)
       else
         c_spin_fock_matrix_mo_oo(m,i) = 0.d0
       endif
-      
-      do e=1,spin_vir_num
-        c_spin_fock_matrix_mo_oo(m,i) = c_spin_fock_matrix_mo_oo(m,i) + 0.5d0*t1_cc(i,e)*spin_fock_matrix_mo_ov(m,e)
-      end do
-      
-      do e=1,spin_vir_num
-        do n=1,spin_occ_num
-          c_spin_fock_matrix_mo_oo(m,i) = c_spin_fock_matrix_mo_oo(m,i) + t1_cc(n,e)*OOOV(m,n,i,e)
-        end do
-      end do
-      
-      do e=1,spin_vir_num
-        do n=1,spin_occ_num
-          do f=1,spin_vir_num
-            c_spin_fock_matrix_mo_oo(m,i) = c_spin_fock_matrix_mo_oo(m,i) + 0.5d0*taus(i,n,e,f)*OOVV(m,n,e,f)
-          end do
-        end do
-      end do
-      
     end do
   end do
   
+      
+  do e=1,spin_vir_num
+    do i=1,spin_occ_num
+      do n=1,spin_occ_num
+        do m=1,spin_occ_num
+          c_spin_fock_matrix_mo_oo(m,i) = c_spin_fock_matrix_mo_oo(m,i) +&
+              OOOV(m,n,i,e) * t1_cc(n,e)
+        end do
+      end do
+    end do
+  end do
+      
+  call dgemm('N','T', &
+    spin_occ_num, spin_occ_num, spin_vir_num, &
+    0.5d0, spin_fock_matrix_mo_ov, spin_occ_num, &
+    t1_cc, spin_occ_num, &
+    1.d0, c_spin_fock_matrix_mo_oo, spin_occ_num)
+
+      
+  call dgemm('N','T',                                                &
+      spin_occ_num, spin_occ_num, spin_vir_num*spin_occ_num*spin_vir_num,&
+      0.5d0, OOVV, spin_occ_num,                                     &
+      taus, spin_occ_num,                                            &
+      1.d0, c_spin_fock_matrix_mo_oo, spin_occ_num)
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, c_spin_fock_matrix_mo_oo_transp, (spin_occ_num,spin_occ_num) ]
@@ -86,7 +91,7 @@ BEGIN_PROVIDER [ double precision, c_spin_fock_matrix_mo_vv, (spin_vir_num,spin_
   
   integer                        :: i,j,m,n
   integer                        :: a,b,e,f
-
+  
   do a=1,spin_vir_num
     do e=1,spin_vir_num
       
@@ -96,25 +101,32 @@ BEGIN_PROVIDER [ double precision, c_spin_fock_matrix_mo_vv, (spin_vir_num,spin_
         c_spin_fock_matrix_mo_vv(a,e) = 0.d0
       endif
       
-      do m=1,spin_occ_num
-        c_spin_fock_matrix_mo_vv(a,e) = c_spin_fock_matrix_mo_vv(a,e) - 0.5d0*t1_cc(m,a)*spin_fock_matrix_mo_ov(m,e)
-      end do
-      
-      do m=1,spin_occ_num
-        do f=1,spin_vir_num
-          c_spin_fock_matrix_mo_vv(a,e) = c_spin_fock_matrix_mo_vv(a,e) + t1_cc(m,f)*OVVV(m,a,f,e)
-        end do
-      end do
-      
-      do m=1,spin_occ_num
-        do n=1,spin_occ_num
-          do f=1,spin_vir_num
-            c_spin_fock_matrix_mo_vv(a,e) = c_spin_fock_matrix_mo_vv(a,e) - 0.5d0*taus(m,n,a,f)*OOVV(m,n,e,f)
-          end do
-        end do
-      end do
-      
     end do
   end do
+  
+  call dgemm('T','N',                                                &
+      spin_vir_num, spin_vir_num, spin_occ_num,                      &
+      -0.5d0, t1_cc, size(t1_cc,1),                                  &
+      spin_fock_matrix_mo_ov, size(spin_fock_matrix_mo_ov,1),        &
+      1.d0, c_spin_fock_matrix_mo_vv, size(c_spin_fock_matrix_mo_vv,1))
+
+  do e=1,spin_vir_num
+    do a=1,spin_vir_num
+      do f=1,spin_vir_num
+        do m=1,spin_occ_num
+          c_spin_fock_matrix_mo_vv(a,e) = c_spin_fock_matrix_mo_vv(a,e) + &
+            t1_cc(m,f)*OVVV(m,a,f,e)
+        end do
+      end do
+    end do
+  end do
+
+
+  call dgemm('T','N',                                                &
+      spin_vir_num, spin_vir_num, spin_occ_num*spin_occ_num*spin_vir_num,&
+      -0.5d0, taus, spin_occ_num*spin_occ_num*spin_vir_num,          &
+      OOVV, spin_occ_num*spin_occ_num*spin_vir_num,                  &
+      1.d0, c_spin_fock_matrix_mo_vv, spin_vir_num)
+  
 END_PROVIDER
 
diff --git a/devel/cc/curly_W.irp.f b/devel/cc/curly_W.irp.f
index 256011a..65b92de 100644
--- a/devel/cc/curly_W.irp.f
+++ b/devel/cc/curly_W.irp.f
@@ -33,7 +33,7 @@ BEGIN_PROVIDER [ double precision, cWoooo, (spin_occ_num,spin_occ_num,spin_occ_n
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, cWovvo_prime, (spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num) ]
+BEGIN_PROVIDER [ double precision, cWovvo, (spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num) ]
  implicit none
  BEGIN_DOC
  ! Curly W in occ-vir-vir-occ block
@@ -43,46 +43,66 @@ BEGIN_PROVIDER [ double precision, cWovvo_prime, (spin_occ_num,spin_vir_num,spin
   double precision               :: x
 
 
-  cWovvo_prime(:,:,:,:) = OVVO_prime(:,:,:,:)
+  cWovvo(:,:,:,:) = OVVO(:,:,:,:)
+
+  call dgemm('N','T', &
+    spin_occ_num*spin_vir_num*spin_vir_num, spin_occ_num, spin_vir_num, &
+    1.d0, OVVV, spin_occ_num*spin_vir_num*spin_vir_num, &
+    t1_cc, spin_occ_num, &
+    1.d0, cWovvo, spin_occ_num*spin_vir_num*spin_vir_num)
 
-  do f=1,spin_vir_num
-    do j=1,spin_occ_num
-      do b=1,spin_vir_num
-        do e=1,spin_vir_num
-          do m=1,spin_occ_num
-            cWovvo_prime(m,e,b,j) = cWovvo_prime(m,e,b,j) + t1_cc(j,f)*OVVV(m,b,e,f)
-          end do
-        end do
-      end do
-    end do
-  end do
 
   do j=1,spin_occ_num
     do e=1,spin_vir_num
-      do b=1,spin_vir_num
-        do n=1,spin_occ_num
-          do m=1,spin_occ_num
-            cWovvo_prime(m,e,b,j) = cWovvo_prime(m,e,b,j) - t1_cc(n,b)*OOVO(m,n,e,j)
-          end do
-        end do
-      end do
+      call dgemm('N','N', &
+        spin_occ_num, spin_vir_num, spin_occ_num, &
+        -1.d0, OOVO(1,1,e,j), spin_occ_num, &
+        t1_cc, spin_occ_num, &
+        1.d0, cWovvo(1,1,e,j), spin_occ_num)
     end do
   end do
 
+  double precision, allocatable :: tmp(:,:,:)
+  allocate ( tmp(spin_occ_num,spin_vir_num,spin_vir_num) )
   do j=1,spin_occ_num
     do f=1,spin_vir_num
       do b=1,spin_vir_num
         do n=1,spin_occ_num
-          x = 0.5d0*t2_cc(j,n,f,b)  + t1_cc(j,f)*t1_cc(n,b)
-          do e=1,spin_vir_num
-            do m=1,spin_occ_num
-              cWovvo_prime(m,e,b,j) = cWovvo_prime(m,e,b,j) - x *OOVV(m,n,e,f)
-            end do
-          end do
+          tmp(n,b,f) = - 0.5d0*t2_cc(j,n,f,b) - t1_cc(j,f)*t1_cc(n,b)
+        enddo
+      enddo
+    enddo
+    do f=1,spin_vir_num
+      do e=1,spin_vir_num
+        call dgemm('N','N', &
+          spin_occ_num, spin_vir_num, spin_occ_num, &
+          1.d0, OOVV(1,1,e,f), spin_occ_num, &
+          tmp(1,1,f), spin_occ_num, &
+          1.d0, cWovvo(1,1,e,j), spin_occ_num)
+      enddo
+    enddo
+  enddo
+  
+END_PROVIDER
+
+  
+BEGIN_PROVIDER [ double precision, cWovvo_prime, (spin_occ_num,spin_vir_num,spin_vir_num,spin_occ_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Curly W in occ-vir-vir-occ block
+ END_DOC
+  integer                        :: m,b,e,j
+
+  do j=1,spin_occ_num
+    do b=1,spin_vir_num
+      do e=1,spin_vir_num
+        do m=1,spin_occ_num
+          cWovvo_prime(m,e,b,j) = cWovvo(m,b,e,j)
         end do
       end do
     end do
   end do
+
   
 END_PROVIDER
 
diff --git a/devel/cc/r1_cc.irp.f b/devel/cc/r1_cc.irp.f
index 48d48cc..683fa6c 100644
--- a/devel/cc/r1_cc.irp.f
+++ b/devel/cc/r1_cc.irp.f
@@ -28,24 +28,16 @@ BEGIN_PROVIDER [ double precision, r1_cc, (spin_occ_num,spin_vir_num) ]
       c_spin_fock_matrix_mo_ov, 1,                                   &
       1.d0, r1_cc, 1)
   
+  double precision, external :: u_dot_v
+
   do f=1,spin_vir_num
     do a=1,spin_vir_num
       do i=1,spin_occ_num
-        do n=1,spin_occ_num
-          r1_cc(i,a) = r1_cc(i,a) - t1_cc(n,f)*OVOV(n,a,i,f)
-        end do
+        r1_cc(i,a) = r1_cc(i,a) - u_dot_v(t1_cc(1,f), OVOV(1,a,i,f),spin_occ_num)
       end do
     end do
-  end do
 
-  call dgemv('T', spin_occ_num*spin_vir_num,                         &
-      spin_occ_num*spin_vir_num,                                     &
-      -1.d0, OVOV, spin_occ_num*spin_vir_num,                        &
-      t1_cc, 1,                                                      &
-      1.d0, r1_cc, 1)
-
-  do e=1,spin_vir_num
-    do f=1,spin_vir_num
+    do e=1,spin_vir_num
       call dgemm('N','N', &
         spin_occ_num, spin_vir_num, spin_occ_num, &
         -0.5d0, t2_cc(1,1,e,f), spin_occ_num, &