Nsryjdtyk

This question is a part of my research. I am trying to multiply two N X N matrices which are stored diagonally. I implemented the following algorithm for this purpose which is working correctly but the problem is it is taking way too much time. For example, To multiply 1000 x 1000 matrices, it is taking almost 26 minutes. I am giving an small example for better understanding of what I am trying to accomplish.

Matrix A = [ 1 2
3 4
],
Matrix B =[ 5 6
7 8],
diagArrayA= [1 4 2 3] (Elements stored diagonally in the order main diagonal - super diagonal -sub diagonal),
diagArrayB = [5 8 6 7] (Same as above),

Result = A X B = [ 19 50 22 43] (Correct output)
The code for multiply function is followed: (The problem seems to be in this function). I will appreciate any comment on the problem/code.

void multiplyDiag(vector <float> diagArrayA, vector <float> diagArrayB){





             cout << "function called " << endl;

           //please declare as array





          // float *a=cutDiagElement(getArray(diagArrayA,-1),0,0);

         //   cout << "Value : " << a[0];



       vector <float> result;

        result.resize(n*n);

      float total_time=0.0;

      // const int sum=n*(n+1)/2;

for(int k=0; k<=n-1;k++)
{ //1

             int start_index=k*n-k*(k-1)/2;



             //cout << start_index;



           //  int stop_index=start_index+n-k-1;



          for(int i=k+1; i<=n-1; i++){ //2



                vector <float> diag_a1=cutDiagElement(getArray(diagArrayA,k-i),0,k);

                // cout << diag_a1[0] <<endl;



                 vector <float> diag_b1=cutDiagElement(getArray(diagArrayB,i),0,0);  



                 int length1=diag_b1.size();

                  //  int length1=n-i;

                    // length1=length1-0-0;

                  //cout << length1;

                   float start1=clock();

                  for(int j=0; j<=length1-1;j++)

                      { //3

                              //cout << "OK";



                     result[start_index+i-k+j]+=diag_a1[j]*diag_b1[j];









                     } //3



                 float end1=clock();

                  total_time=total_time+(end1-start1);







                  vector <float> diag_a2=cutDiagElement(getArray(diagArrayA, i), 0, 0);



                  vector <float> diag_b2= cutDiagElement(getArray(diagArrayB,k-i), k, 0);

                   //  cout << diag_b2[0];

                //  int length2=diag_a2.size();



                   float start2=clock();

                   for(int j=0; j<=length1-1;j++)

                       { //4





                              result[start_index+j]+=diag_a2[j]*diag_b2[j];

                               //  cout<< multiply_result[start_index+j] << endl;



                        }    //4 



                 float end2=clock();



                total_time=total_time+(end2-start2);





       } //2



   for(int i=0; i<=k; i++)

     { //5

             vector <float> diag_a3=cutDiagElement(getArray(diagArrayA, i), 0, k-i);



              vector <float> diag_b3 = cutDiagElement(getArray(diagArrayB, k-i), i, 0);



              int length3=diag_a3.size();

             //   int length3=n-i;

               //  length3=length3-(k-i);



              float start3=clock();



                for(int j=0; j<=length3-1;j++){ //6



                                           result[start_index+j]+=diag_a3[j]*diag_b3[j];



                                             //cout<< multiply_result[start_index+j] << endl;



                              }  //6



           float end3=clock();

          total_time=total_time+(end3-start3);





    } //5

}//1

// computing sub diagonal

for(int k=1; k<=n-1; k++)
{ //7

    int start_index=n*(n+1)/2+(k-1)*n-k*(k-1)/2;

      // int stop_index=start_index+n-k-1;

    for(int i=k+1; i<=n-1; i++)

     { //8

              vector <float> diag_a1=cutDiagElement(getArray(diagArrayA, -i), 0, 0);

               vector <float> diag_b1= cutDiagElement(getArray(diagArrayB, i-k), 0, k);

               int length1=diag_b1.size();

              // int length1=n-(i-k);

                //   length1=length1-k;

               float start4=clock();

               for(int j=0; j<=length1-1;j++){ //9



                          result[start_index+i+j-k]+=diag_a1[j]*diag_b1[j];





                    } //9



                float end4=clock();

                  total_time=total_time+(end4-start4);



                 vector <float> diag_a2=cutDiagElement(getArray(diagArrayA, i-k), k, 0);



                 vector <float> diag_b2= cutDiagElement(getArray(diagArrayB,-i), 0, 0);

                // int length2=diag_a2.size();



                float start5=clock();



                 for(int j=0; j<=length1-1;j++)

                           { //10



                             result[start_index+j]+=diag_a2[j]*diag_b2[j];



                          }//10  



                   float end5=clock();



               total_time=total_time+(end5-start5);





     } //8



   for(int i=0; i<=k; i++)

               { //11

                    vector <float> diag_a3=cutDiagElement(getArray(diagArrayA, -i), k-i, 0);



                     vector <float> diag_b3 = cutDiagElement(getArray(diagArrayB, i-k), 0, i);



                      int length3=diag_a3.size();

                     //     int length3=n-i;

                       //     length3=length3-(k-i);



                    float start6=clock();



                    for(int j=0; j<=length3-1;j++){ //12



                                           result[start_index+j]+=diag_a3[j]*diag_b3[j];



                              }  //12

                   float end6=clock();

                  total_time=total_time+(end6-start6);





                } //11

} //7

cout << "time :" << (total_time) / CLOCKS_PER_SEC*1000 <

           cout << result[x] << endl;



           }

}