Changeset b055de9 in git


Ignore:
Timestamp:
May 9, 2005, 11:39:28 AM (18 years ago)
Author:
Michael Brickenstein <bricken@…>
Branches:
(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
Children:
81306ac9be57a712101e2e5b63ade473de30abd8
Parents:
f4adfcbbf4e6a189938d7ffb03fbf5d42b6a0250
Message:
*bricken: gauss externalized


git-svn-id: file:///usr/local/Singular/svn/trunk@8101 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
kernel
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • kernel/Makefile.in

    rf4adfcb rb055de9  
    113113    pDebug.cc pInline2.cc pInline1.cc pInline0.cc \
    114114    pShallowCopyDelete.cc fast_mult.cc digitech.cc\
    115     tgb.cc
     115    tgb.cc tgbgauss.cc
    116116
    117117# normal C source files

  • kernel/tgb.cc

    rf4adfcb rb055de9  
    1212#include "tgb.h"
    1313#include "tgb_internal.h"
     14#include "tgbgauss.h"
    1415static const int bundle_size=100;
    1516#if 1
    … …  
    23222323  assume(F_minus==NULL);
    23232324}
    2324 tgb_matrix::tgb_matrix(int i, int j){
    2325   n=(number**) omalloc(i*sizeof (number*));;
    2326   int z;
    2327   int z2;
    2328   for(z=0;z<i;z++)
    2329   {
    2330     n[z]=(number*)omalloc(j*sizeof(number));
    2331     for(z2=0;z2<j;z2++)
    2332     {
    2333       n[z][z2]=nInit(0);
    2334     }
    2335   }
    2336   this->columns=j;
    2337   this->rows=i;
    2338   free_numbers=FALSE;
    2339 }
    2340 tgb_matrix::~tgb_matrix(){
    2341   int z;
    2342   for(z=0;z<rows;z++)
    2343   {
    2344     if(n[z])
    2345     {
    2346       if(free_numbers)
    2347       {
    2348         int z2;
    2349         for(z2=0;z2<columns;z2++)
    2350         {
    2351           nDelete(&(n[z][z2]));
    2352         }
    2353       }
    2354       omfree(n[z]);
    2355     }
    2356   }
    2357   omfree(n);
    2358 }
    2359 void tgb_matrix::print(){
    2360   int i;
    2361   int j;
    2362   Print("\n");
    2363   for(i=0;i<rows;i++)
    2364   {
    2365     Print("(");
    2366     for(j=0;j<columns;j++)
    2367     {
    2368       StringSetS("");
    2369       n_Write(n[i][j],currRing);
    2370       Print(StringAppendS(""));
    2371       Print("\t");
    2372     }
    2373     Print(")\n");
    2374   }
    2375 }
    2376 //transfers ownership of n to the matrix
    2377 void tgb_matrix::set(int i, int j, number n){
    2378   assume(i<rows);
    2379   assume(j<columns);
    2380   this->n[i][j]=n;
    2381 }
    2382 int tgb_matrix::get_rows(){
    2383   return rows;
    2384 }
    2385 int tgb_matrix::get_columns(){
    2386   return columns;
    2387 }
    2388 number tgb_matrix::get(int i, int j){
    2389   assume(i<rows);
    2390   assume(j<columns);
    2391   return n[i][j];
    2392 }
    2393 BOOLEAN tgb_matrix::is_zero_entry(int i, int j){
    2394   return (nIsZero(n[i][j]));
    2395 }
    2396 void tgb_matrix::perm_rows(int i, int j){
    2397   number* h;
    2398   h=n[i];
    2399   n[i]=n[j];
    2400   n[j]=h;
    2401 }
    2402 int tgb_matrix::min_col_not_zero_in_row(int row){
    2403   int i;
    2404   for(i=0;i<columns;i++)
    2405   {
    2406     if(!(nIsZero(n[row][i])))
    2407       return i;
    2408   }
    2409   return columns;//error code
    2410 }
    2411 int tgb_matrix::next_col_not_zero(int row,int pre){
    2412   int i;
    2413   for(i=pre+1;i<columns;i++)
    2414   {
    2415     if(!(nIsZero(n[row][i])))
    2416       return i;
    2417   }
    2418   return columns;//error code
    2419 }
    2420 BOOLEAN tgb_matrix::zero_row(int row){
    2421   int i;
    2422   for(i=0;i<columns;i++)
    2423   {
    2424     if(!(nIsZero(n[row][i])))
    2425       return FALSE;
    2426   }
    2427   return TRUE;
    2428 }
    2429 int tgb_matrix::non_zero_entries(int row){
    2430   int i;
    2431   int z=0;
    2432   for(i=0;i<columns;i++)
    2433   {
    2434     if(!(nIsZero(n[row][i])))
    2435       z++;
    2436   }
    2437   return z;
    2438 }
    2439 //row add_to=row add_to +row summand*factor
    2440 void tgb_matrix::add_lambda_times_row(int add_to,int summand,number factor){
    2441   int i;
    2442   for(i=0;i<columns;i++){
    2443     if(!(nIsZero(n[summand][i])))
    2444     {
    2445       number n1=n[add_to][i];
    2446       number n2=nMult(factor,n[summand][i]);
    2447       n[add_to][i]=nAdd(n1,n2);
    2448       nDelete(&n1);
    2449       nDelete(&n2);
    2450     }
    2451   }
    2452 }
    2453 void tgb_matrix::mult_row(int row,number factor){
    2454   if (nIsOne(factor))
    2455     return;
    2456   int i;
    2457   for(i=0;i<columns;i++){
    2458     if(!(nIsZero(n[row][i])))
    2459     {
    2460       number n1=n[row][i];
    2461       n[row][i]=nMult(n1,factor);
    2462       nDelete(&n1);
    2463     }
    2464   }
    2465 }
    2466 void tgb_matrix::free_row(int row, BOOLEAN free_non_zeros){
    2467   int i;
    2468   for(i=0;i<columns;i++)
    2469     if((free_non_zeros)||(!(nIsZero(n[row][i]))))
    2470       nDelete(&(n[row][i]));
    2471   omfree(n[row]);
    2472   n[row]=NULL;
    2473 }
    2474 
    2475 
    2476 tgb_sparse_matrix::tgb_sparse_matrix(int i, int j, ring rarg){
    2477   mp=(mac_poly*) omalloc(i*sizeof (mac_poly));;
    2478   int z;
    2479   int z2;
    2480   for(z=0;z<i;z++)
    2481   {
    2482     mp[z]=NULL;
    2483   }
    2484   this->columns=j;
    2485   this->rows=i;
    2486   free_numbers=FALSE;
    2487   r=rarg;
    2488 }
    2489 tgb_sparse_matrix::~tgb_sparse_matrix(){
    2490   int z;
    2491   for(z=0;z<rows;z++)
    2492   {
    2493     if(mp[z])
    2494     {
    2495       if(free_numbers)
    2496       {
    2497         mac_destroy(mp[z]);
    2498       }
    2499       else {
    2500         while(mp[z])
    2501         {
    2502          
    2503           mac_poly next=mp[z]->next;
    2504           delete mp[z];
    2505           mp[z]=next;
    2506         }
    2507       }
    2508     }
    2509   }
    2510   omfree(mp);
    2511 }
    2512 int row_cmp_gen(const void* a, const void* b){
    2513   const mac_poly ap= *((mac_poly*) a);
    2514   const mac_poly bp=*((mac_poly*) b);
    2515   if (ap==NULL) return 1;
    2516   if (bp==NULL) return -1;
    2517   if (ap->exp<bp->exp) return -1;
    2518   return 1;
    2519 }
    2520 void tgb_sparse_matrix::sort_rows(){
    2521   qsort(mp,rows,sizeof(mac_poly),row_cmp_gen);
    2522 }
    2523 void tgb_sparse_matrix::print(){
    2524   int i;
    2525   int j;
    2526   Print("\n");
    2527   for(i=0;i<rows;i++)
    2528   {
    2529     Print("(");
    2530     for(j=0;j<columns;j++)
    2531     {
    2532       StringSetS("");
    2533       number n=get(i,j);
    2534       n_Write(n,currRing);
    2535       Print(StringAppendS(""));
    2536       Print("\t");
    2537     }
    2538     Print(")\n");
    2539   }
    2540 }
    2541 //transfers ownership of n to the matrix
    2542 void tgb_sparse_matrix::set(int i, int j, number n){
    2543   assume(i<rows);
    2544   assume(j<columns);
    2545   mac_poly* set_this=&mp[i];
    2546   //  while(((*set_this)!=NULL)&&((*set_this)­>exp<j))
    2547   while(((*set_this)!=NULL) && ((*set_this)->exp<j))
    2548     set_this=&((*set_this)->next);
    2549 
    2550   if (((*set_this)==NULL)||((*set_this)->exp>j))
    2551   {
    2552     if (nIsZero(n)) return;
    2553     mac_poly old=(*set_this);
    2554     (*set_this)=new mac_poly_r();
    2555     (*set_this)->exp=j;
    2556     (*set_this)->coef=n;
    2557     (*set_this)->next=old;
    2558     return;
    2559   }
    2560   assume((*set_this)->exp==j);
    2561   if(!nIsZero(n))
    2562   {
    2563     nDelete(&(*set_this)->coef);
    2564     (*set_this)->coef=n;
    2565   }
    2566   else
    2567   {
    2568     nDelete(&(*set_this)->coef);
    2569     mac_poly dt=(*set_this);
    2570     (*set_this)=dt->next;
    2571     delete dt;
    2572    
    2573    
    2574   }
    2575   return;
    2576 }
    2577 
    2578 
    2579 
    2580 int tgb_sparse_matrix::get_rows(){
    2581   return rows;
    2582 }
    2583 int tgb_sparse_matrix::get_columns(){
    2584   return columns;
    2585 }
    2586 number tgb_sparse_matrix::get(int i, int j){
    2587   assume(i<rows);
    2588   assume(j<columns);
    2589   mac_poly r=mp[i];
    2590   while((r!=NULL)&&(r->exp<j))
    2591     r=r->next;
    2592   if ((r==NULL)||(r->exp>j))
    2593   {
    2594     number n=nInit(0);
    2595     return n;
    2596   }
    2597   assume(r->exp==j);
    2598   return r->coef;
    2599 }
    2600 BOOLEAN tgb_sparse_matrix::is_zero_entry(int i, int j){
    2601   assume(i<rows);
    2602   assume(j<columns);
    2603   mac_poly r=mp[i];
    2604   while((r!=NULL)&&(r->exp<j))
    2605     r=r->next;
    2606   if ((r==NULL)||(r->exp>j))
    2607   {
    2608     return TRUE;
    2609   }
    2610   assume(!nIsZero(r->coef));
    2611   assume(r->exp==j);
    2612   return FALSE;
    2613  
    2614 }
    2615 
    2616 int tgb_sparse_matrix::min_col_not_zero_in_row(int row){
    2617   if(mp[row]!=NULL)
    2618   {
    2619     assume(!nIsZero(mp[row]->coef));
    2620     return mp[row]->exp;
    2621   }
    2622   else
    2623 
    2624  
    2625   return columns;//error code
    2626 }
    2627 int tgb_sparse_matrix::next_col_not_zero(int row,int pre){ 
    2628   mac_poly r=mp[row];
    2629   while((r!=NULL)&&(r->exp<=pre))
    2630     r=r->next;
    2631   if(r!=NULL)
    2632   {
    2633     assume(!nIsZero(r->coef));
    2634     return r->exp;
    2635   }
    2636   return columns;//error code
    2637 }
    2638 BOOLEAN tgb_sparse_matrix::zero_row(int row){
    2639   assume((mp[row]==NULL)||(!nIsZero(mp[row]->coef)));
    2640   if (mp[row]==NULL)
    2641     return TRUE;
    2642   else
    2643     return FALSE;
    2644 }
    2645 void tgb_sparse_matrix::row_normalize(int row){
    2646   if (!rField_has_simple_inverse(r))  /* Z/p, GF(p,n), R, long R/C */
    2647   {
    2648     mac_poly m=mp[row];
    2649         while (m!=NULL)
    2650         {
    2651           if (currRing==r) {nTest(m->coef);}
    2652           n_Normalize(m->coef,r);
    2653           m=m->next;
    2654         }
    2655   }
    2656 }
    2657 void tgb_sparse_matrix::row_content(int row){
    2658  
    2659   mac_poly ph=mp[row];
    2660   number h,d;
    2661   mac_poly p;
    2662 
    2663   if(TEST_OPT_CONTENTSB) return;
    2664   if(ph->next==NULL)
    2665   {
    2666     nDelete(&ph->coef);
    2667     ph->coef=nInit(1);
    2668   }
    2669   else
    2670   {
    2671     nNormalize(ph->coef);
    2672     if(!nGreaterZero(ph->coef)) {
    2673       //ph = pNeg(ph);
    2674       p=ph;
    2675       while(p)
    2676       {
    2677         p->coef=nNeg(p->coef);
    2678         p=p->next;
    2679       }
    2680     }
    2681    
    2682     h=nCopy(ph->coef);
    2683     p = ph->next;
    2684    
    2685     while (p!=NULL)
    2686     {
    2687       nNormalize(p->coef);
    2688       d=nGcd(h,p->coef,currRing);
    2689       nDelete(&h);
    2690       h = d;
    2691       if(nIsOne(h))
    2692       {
    2693         break;
    2694       }
    2695       p=p->next;
    2696     }
    2697     p = ph;
    2698     //number tmp;
    2699     if(!nIsOne(h))
    2700     {
    2701       while (p!=NULL)
    2702       {
    2703    
    2704         d = nIntDiv(pGetCoeff(p),h);
    2705         nDelete(&p->coef);
    2706         p->coef=d;
    2707         p=p->next;
    2708       }
    2709     }
    2710     nDelete(&h);
    2711 
    2712   }
    2713 }
    2714 int tgb_sparse_matrix::non_zero_entries(int row){
    2715 
    2716   return mac_length(mp[row]);
    2717 }
    2718 //row add_to=row add_to +row summand*factor
    2719 void tgb_sparse_matrix::add_lambda_times_row(int add_to,int summand,number factor){
    2720   mp[add_to]= mac_p_add_ff_qq(mp[add_to], factor,mp[summand]);
    2721 
    2722 }
    2723 void tgb_sparse_matrix::mult_row(int row,number factor){
    2724   if (nIsZero(factor))
    2725   {
    2726     mac_destroy(mp[row]);
    2727     mp[row]=NULL;
    2728    
    2729     return;
    2730   }
    2731   if(nIsOne(factor))
    2732     return;
    2733   mac_mult_cons(mp[row],factor);
    2734 }
    2735 void tgb_sparse_matrix::free_row(int row, BOOLEAN free_non_zeros){
    2736   if(free_non_zeros)
    2737     mac_destroy(mp[row]);
    2738   else
    2739   {
    2740     while(mp[row])
    2741     {
    2742      
    2743       mac_poly next=mp[row]->next;
    2744       delete mp[row];
    2745       mp[row]=next;
    2746     }
    2747   }
    2748   mp[row]=NULL;
    2749 }
     2325
    27502326//transfers ownership of m to mat
    27512327void init_with_mac_poly(tgb_sparse_matrix* mat, int row, mac_poly m){
    … …  
    27792355}
    27802356
    2781 void simple_gauss(tgb_sparse_matrix* mat, calc_dat* c){
    2782   int col, row;
    2783   int* row_cache=(int*) omalloc(mat->get_rows()*sizeof(int));
    2784   col=0;
    2785   row=0;
    2786   int i;
    2787   int pn=mat->get_rows();
    2788   int matcol=mat->get_columns();
    2789   int* area=(int*) omalloc(sizeof(int)*((matcol-1)/bundle_size+1));
    2790   const int max_area_index=(matcol-1)/bundle_size;
    2791     //rows are divided in areas
    2792   //if row begins with columns col, it is located in [area[col/bundle_size],area[col/bundle_size+1]-1]
    2793   assume(pn>0);
    2794   //first clear zeroes
    2795   for(i=0;i<pn;i++)
    2796   {
    2797     if(mat->zero_row(i))
    2798     {
    2799       mat->perm_rows(i,pn-1);
    2800       pn--;
    2801       if(i!=pn){i--;}
    2802     }
    2803  
    2804   }
    2805   mat->sort_rows();
    2806   for(i=0;i<pn;i++)
    2807   {
    2808       row_cache[i]=mat->min_col_not_zero_in_row(i);
    2809       // Print("row_cache:%d\n",row_cache[i]);
    2810   }
    2811   int last_area=-1;
    2812   for(i=0;i<pn;i++)
    2813   {
    2814     int this_area=row_cache[i]/bundle_size;
    2815     assume(this_area>=last_area);
    2816     if(this_area>last_area)
    2817     {
    2818       int j;
    2819       for(j=last_area+1;j<=this_area;j++)
    2820         area[j]=i;
    2821       last_area=this_area;
    2822     }
    2823   }
    2824   for(i=last_area+1;i<=max_area_index;i++)
    2825   {
    2826     area[i]=pn;
    2827   }
    2828   while(row<pn-1){
    2829     //row is the row where pivot should be
    2830     // row== pn-1 means we have only to act on one row so no red nec.
    2831     //we assume further all rows till the pn-1 row are non-zero
    2832    
    2833     //select column
    2834    
    2835     //col=mat->min_col_not_zero_in_row(row);
    2836     int max_in_area;
    2837     {
    2838       int tai=row_cache[row]/bundle_size;
    2839       assume(tai<=max_area_index);
    2840       if(tai==max_area_index)
    2841         max_in_area=pn-1;
    2842       else
    2843         max_in_area=area[tai+1]-1;
    2844     }
    2845     assume(row_cache[row]==mat->min_col_not_zero_in_row(row));
    2846     col=row_cache[row];
    2847 
    2848    
    2849     assume(col!=matcol);
    2850     int found_in_row;
    2851    
    2852     found_in_row=row;
    2853     BOOLEAN must_reduce=FALSE;
    2854     assume(pn<=mat->get_rows());
    2855     for(i=row+1;i<=max_in_area;i++){
    2856       int first;//=mat->min_col_not_zero_in_row(i);
    2857       assume(row_cache[i]==mat->min_col_not_zero_in_row(i));
    2858       first=row_cache[i];
    2859       assume(first!=matcol);
    2860       if(first<col)
    2861       {
    2862         col=first;
    2863         found_in_row=i;
    2864         must_reduce=FALSE;
    2865       }
    2866       else {
    2867         if(first==col)
    2868           must_reduce=TRUE;
    2869       }
    2870     }
    2871     //select pivot
    2872     int act_l=nSize(mat->get(found_in_row,col))*mat->non_zero_entries(found_in_row);
    2873     if(must_reduce)
    2874     {
    2875       for(i=found_in_row+1;i<=max_in_area;i++){
    2876         assume(mat->min_col_not_zero_in_row(i)>=col);
    2877         int first;
    2878         assume(row_cache[i]==mat->min_col_not_zero_in_row(i));
    2879         first=row_cache[i];
    2880         assume(first!=matcol);
    2881         //      if((!(mat->is_zero_entry(i,col)))&&(mat->non_zero_entries(i)<act_l))
    2882         int nz;
    2883         if((row_cache[i]==col)&&((nz=nSize(mat->get(i,col))*mat->non_zero_entries(i))<act_l))
    2884         {
    2885           found_in_row=i;
    2886         act_l=nz;
    2887         }
    2888        
    2889       }
    2890     }
    2891     mat->perm_rows(row,found_in_row);
    2892     int h=row_cache[row];
    2893     row_cache[row]=row_cache[found_in_row];
    2894     row_cache[found_in_row]=h;
    2895 
    2896 
    2897 
    2898     if(!must_reduce){
    2899       row++;
    2900       continue;
    2901     }
    2902     //reduction
    2903     //must extract content and normalize here
    2904     mat->row_content(row);
    2905     mat->row_normalize(row);
    2906 
    2907     //for(i=row+1;i<pn;i++){
    2908     for(i=max_in_area;i>row;i--)
    2909     {
    2910       int col_area_index=col/bundle_size;
    2911       assume(col_area_index<=max_area_index);
    2912       assume(mat->min_col_not_zero_in_row(i)>=col);
    2913       int first;
    2914       assume(row_cache[i]==mat->min_col_not_zero_in_row(i));
    2915       first=row_cache[i];
    2916       assume(first!=matcol);
    2917       if(row_cache[i]==col)
    2918       {
    2919        
    2920         number c1=mat->get(i,col);
    2921         number c2=mat->get(row,col);
    2922         number n1=c1;
    2923         number n2=c2;
    2924 
    2925         ksCheckCoeff(&n1,&n2);
    2926         //nDelete(&c1);
    2927         n1=nNeg(n1);
    2928         mat->mult_row(i,n2);
    2929         mat->add_lambda_times_row(i,row,n1);
    2930         nDelete(&n1);
    2931         nDelete(&n2);
    2932         assume(mat->is_zero_entry(i,col));
    2933         row_cache[i]=mat->min_col_not_zero_in_row(i);
    2934         assume(mat->min_col_not_zero_in_row(i)>col);
    2935         if(row_cache[i]==matcol)
    2936         {
    2937           int index;
    2938           index=i;
    2939           int last_in_area;
    2940           int this_cai=col_area_index;
    2941           while(this_cai<max_area_index)
    2942           {
    2943             last_in_area=area[this_cai+1]-1;
    2944             int h_c=row_cache[last_in_area];
    2945             row_cache[last_in_area]=row_cache[index];
    2946             row_cache[index]=h_c;
    2947             mat->perm_rows(index,last_in_area);
    2948             index=last_in_area;
    2949             this_cai++;
    2950             area[this_cai]--;
    2951           }
    2952           mat->perm_rows(index,pn-1);
    2953           row_cache[index]=row_cache[pn-1];
    2954           row_cache[pn-1]=matcol;
    2955           pn--;
    2956         }
    2957         else
    2958         {
    2959           int index;
    2960           index=i;
    2961           int last_in_area;
    2962           int this_cai=col_area_index;
    2963           int final_cai=row_cache[index]/bundle_size;
    2964           assume(final_cai<=max_area_index);
    2965           while(this_cai<final_cai)
    2966           {
    2967             last_in_area=area[this_cai+1]-1;
    2968             int h_c=row_cache[last_in_area];
    2969             row_cache[last_in_area]=row_cache[index];
    2970             row_cache[index]=h_c;
    2971             mat->perm_rows(index,last_in_area);
    2972             index=last_in_area;
    2973             this_cai++;
    2974             area[this_cai]--;
    2975           }
    2976         }
    2977       }
    2978       else
    2979         assume(mat->min_col_not_zero_in_row(i)>col);
    2980     }
    2981 
    2982 //     for(i=row+1;i<pn;i++)
    2983 //     {
    2984 //       assume(mat->min_col_not_zero_in_row(i)==row_cache[i]);
    2985 //       // if(mat->zero_row(i))
    2986 //       assume(matcol==mat->get_columns());
    2987 //       if(row_cache[i]==matcol)
    2988 //       {
    2989 //      assume(mat->zero_row(i));
    2990 //      mat->perm_rows(i,pn-1);
    2991 //      row_cache[i]=row_cache[pn-1];
    2992 //      row_cache[pn-1]=matcol;
    2993 //      pn--;
    2994 //      if(i!=pn){i--;}
    2995 //       }
    2996 //     }
    2997 #ifdef TGB_DEBUG
    2998   {
    2999   int last=-1;
    3000   for(i=0;i<pn;i++)
    3001   {
    3002     int act=mat->min_col_not_zero_in_row(i);
    3003     assume(act>last);
    3004    
    3005   }
    3006   for(i=pn;i<mat->get_rows();i++)
    3007   {
    3008     assume(mat->zero_row(i));
    3009    
    3010   }
    3011  
    3012 
    3013   }
    3014 #endif
    3015     row++;
    3016   }
    3017   omfree(area);
    3018   omfree(row_cache);
    3019 }
    3020 void simple_gauss2(tgb_matrix* mat){
    3021   int col, row;
    3022   col=0;
    3023   row=0;
    3024   int i;
    3025   int pn=mat->get_rows();
    3026   assume(pn>0);
    3027   //first clear zeroes
    3028 //   for(i=0;i<pn;i++)
    3029 //   {
    3030 //     if(mat->zero_row(i))
    3031 //     {
    3032 //       mat->perm_rows(i,pn-1);
    3033 //       pn--;
    3034 //       if(i!=pn){i--;}
    3035 //     }
    3036 //   }
    3037   while((row<pn-1)&&(col<mat->get_columns())){
    3038     //row is the row where pivot should be
    3039     // row== pn-1 means we have only to act on one row so no red nec.
    3040     //we assume further all rows till the pn-1 row are non-zero
    3041    
    3042     //select column
    3043    
    3044     //    col=mat->min_col_not_zero_in_row(row);
    3045     assume(col!=mat->get_columns());
    3046     int found_in_row=-1;
    3047    
    3048     //    found_in_row=row;
    3049     assume(pn<=mat->get_rows());
    3050     for(i=row;i<pn;i++)
    3051     {
    3052       //    int first=mat->min_col_not_zero_in_row(i);
    3053       //  if(first<col)
    3054       if(!(mat->is_zero_entry(i,col))){
    3055        
    3056         found_in_row=i;
    3057         break;
    3058       }
    3059     }
    3060     if(found_in_row!=-1)
    3061     {
    3062     //select pivot
    3063       int act_l=mat->non_zero_entries(found_in_row);
    3064       for(i=i+1;i<pn;i++)
    3065       {
    3066         int vgl;
    3067         assume(mat->min_col_not_zero_in_row(i)>=col);
    3068         if((!(mat->is_zero_entry(i,col)))&&((vgl=mat->non_zero_entries(i))<act_l))
    3069         {
    3070           found_in_row=i;
    3071           act_l=vgl;
    3072         }
    3073        
    3074       }
    3075       mat->perm_rows(row,found_in_row);
    3076      
    3077      
    3078       //reduction
    3079       for(i=row+1;i<pn;i++){
    3080         assume(mat->min_col_not_zero_in_row(i)>=col);
    3081         if(!(mat->is_zero_entry(i,col)))
    3082         {
    3083          
    3084           number c1=nNeg(nCopy(mat->get(i,col)));
    3085           number c2=mat->get(row,col);
    3086           number n1=c1;
    3087           number n2=c2;
    3088          
    3089           ksCheckCoeff(&n1,&n2);
    3090           nDelete(&c1);
    3091           mat->mult_row(i,n2);
    3092           mat->add_lambda_times_row(i,row,n1);
    3093           assume(mat->is_zero_entry(i,col));
    3094          
    3095         }
    3096         assume(mat->min_col_not_zero_in_row(i)>col);
    3097       }
    3098       row++;
    3099     }
    3100     col++;
    3101     // for(i=row+1;i<pn;i++)
    3102 //     {
    3103 //       if(mat->zero_row(i))
    3104 //       {
    3105 //      mat->perm_rows(i,pn-1);
    3106 //      pn--;
    3107 //      if(i!=pn){i--;}
    3108 //       }
    3109 //     }
    3110 
    3111   }
    3112 }
     2357
    31132358
    31142359static tgb_matrix* build_matrix(poly* p,int p_index,poly* done, int done_index, calc_dat* c){

  • kernel/tgb_internal.h

    rf4adfcb rb055de9  
    3333#define AC_NEW_MIN 2
    3434#define AC_FLATTEN 1
     35
    3536//#define FIND_DETERMINISTIC
    3637//#define REDTAIL_PROT
    … …  
    101102  mp_array_list* next;
    102103};
    103 class tgb_matrix{
    104  private:
    105   number** n;
    106   int columns;
    107   int rows;
    108   BOOLEAN free_numbers;
    109  public:
    110   tgb_matrix(int i, int j);
    111   ~tgb_matrix();
    112   int get_rows();
    113   int get_columns();
    114   void print();
    115   void perm_rows(int i, int j);
    116   void set(int i, int j, number n);
    117   number get(int i, int j);
    118   BOOLEAN is_zero_entry(int i, int j);
    119   void free_row(int row, BOOLEAN free_non_zeros=TRUE);
    120   int min_col_not_zero_in_row(int row);
    121   int next_col_not_zero(int row,int pre);
    122   BOOLEAN zero_row(int row);
    123   void mult_row(int row,number factor);
    124   void add_lambda_times_row(int add_to,int summand,number factor);
    125   int non_zero_entries(int row);
    126 };
     104
    127105class mac_poly_r{
    128106public:
    … …  
    136114
    137115typedef mac_poly_r* mac_poly;
    138 class tgb_sparse_matrix{
    139  private:
    140   ring r;
    141   mac_poly* mp;
    142   int columns;
    143   int rows;
    144   BOOLEAN free_numbers;
    145  public:
    146   void sort_rows();
    147   friend poly free_row_to_poly(tgb_sparse_matrix* mat, int row, poly* monoms, int monom_index);
    148   friend void init_with_mac_poly(tgb_sparse_matrix* mat, int row, mac_poly m);
    149   tgb_sparse_matrix(int i, int j, ring rarg);
    150   ~tgb_sparse_matrix();
    151   int get_rows();
    152   int get_columns();
    153   void print();
    154   void row_normalize(int row);
    155   void row_content(int row);
    156   //  void perm_rows(int i, int j);
    157   void perm_rows(int i, int j){
    158   mac_poly h;
    159   h=mp[i];
    160   mp[i]=mp[j];
    161   mp[j]=h;
    162   }
    163   void set(int i, int j, number n);
    164   number get(int i, int j);
    165   BOOLEAN is_zero_entry(int i, int j);
    166   void free_row(int row, BOOLEAN free_non_zeros=TRUE);
    167   int min_col_not_zero_in_row(int row);
    168   int next_col_not_zero(int row,int pre);
    169   BOOLEAN zero_row(int row);
    170   void mult_row(int row,number factor);
    171   void add_lambda_times_row(int add_to,int summand,number factor);
    172   int non_zero_entries(int row);
    173 };
     116
    174117struct poly_array_list{
    175118  poly* p;
    … …  
    335278static void multi_reduce_step(find_erg & erg, red_object* r, calc_dat* c);
    336279static void finalize_reduction_step(reduction_step* r);
     280
     281
     282
     283
     284
     285
     286
     287
     288mac_poly mac_p_add_ff_qq(mac_poly a, number f,mac_poly b);
     289
     290void mac_mult_cons(mac_poly p,number c);
     291int mac_length(mac_poly p);
     292
     293//contrary to delete on the mac_poly_r, the coefficients are also destroyed here
     294void mac_destroy(mac_poly p);
     295
     296
     297
    337298#endif
Note: See TracChangeset for help on using the changeset viewer.